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Abstract 

This  thesis  describes  the  development  of  a  Fortran  computer  code  which  models  the  interaction  between 
an  incompressible,  potential  flow  and  a  homogeneous,  elastic  structure.  The  boundary  element  technique  was 
chosen  over  finite  differences  and  finite  elements  because  of  its  ability  to  numerically  approximate  both  the  fluid 
and  structural  behavior  with  a  common  definition  of  the  fluid/structure  boundary. 

The  ability  to  accurately  model  solid  and  fluid  boundaries  can  be  quite  important  in  the  fields  of 
aeroelasticity  and  structural  analysis.  The  nature  of  these  boundaries  is  what  determines  the  final  solution  to 
a  problem  of  fluid  flow  past  an  elastic  body.  Of‘?n  the  complexity  of  defining  and  tracking  the  boundary  and 
its  associated  boundary  conditions  has  led  the  user  to  assumptions  of  rigid  bodies,  and  therefore  rigid 
boundaries.  Certainly  the  tasks  of  defining  the  domain  grids  for  finite  difference  and  finite  element  techniques 
have  not  simplified  the  process. 

In  the  computer  code  developed  for  this  thesis  the  fluid  and  structural  governing  equations  are  simulta¬ 
neously  solved  to  determine  the  pressure  about  the  structure  and  the  corresponding  elastic  deformations.  The 
deformations  are  applied  to  the  original  boundary,  resulting  in  a  new  geometry.  This  new  geometry  is  used  to 
recalculate  the  pressure  field  about  the  structure,  and  the  process  is  iterated  until  a  final  steady-state  solution 
is  obtained. 

While  simplifying  assumptions  have  been  made  in  the  execution  of  this  thesis,  the  general  boundary 
element  techique  has  the  ability  to  model  complex,  higher  order  problems.  The  initial  results  obtained  during 
the  course  of  this  work  show  promise,  and  follow-on  studies  are  recommended. 


) 

SOLUTION  OF  POTENTIAL  FLOW  PAST  AN  ELASTIC  BODY 
USING  THE  BOUNDARY  ELEMENT  TECHNIQUE 

I.  Introduction 

This  thesis  investigates  the  simultaneous  solution  of  the  interaction  between  an  incompressible, 
potential  flow  and  a  homogeneous,  elastic  structure  using  the  boundary  element  technique.  The  ability  to 
accurately  model  the  fluid  and  structural  boundaries  can  be  quite  important  in  the  Fields  of  aeroelasticity  as 
the  emphasis  in  aeroelasticity  is; 

.  .  .  phenomena  which  exhibit  appreciable  reciprocal  interaction  (static  or  dynamic)  between 
aerodynamic  forces  and  the  deformations  induced  thereby  in  the  structure  of  a  flying  vehicle,  its 
control  mechanisms,  or  its  propulsion  system.  (1:1) 

For  example,  the  twist  angle-of-attack  of  a  wing  is  composed  of  a  series  of  incremental  twists  each  caused 
by  the  aerodynamic  loading  of  the  previous  geometry  (1:  Eq  6-7).  The  original  wing  geometry  causes 
aerodynamic  loading  which  in  turn  produces  a  twist,  or  deformation,  to  the  wing.  This  deformation  causes 
a  slightly  changed  » jading  which  in  turn  causes  additional  twist,  and  so  on.  The  ability  to  compute  this 
incremental  deformation  determines  the  final  solution  to  a  problem  of  fluid  flow  past  an  elastic  body.  The 
traditional  approach  to  aeroelasticity  is  to  perform  a  modal  analysis  of  the  structure  assuming  structural 
shapes  and  determining  the  aerodynamic  loads  due  to  those  shapes.  The  approach  taken  in  this  research 
presents  a  very  different  one  in  that  the  fluid  and  structure  governing  equations  are  simultaneously  solved 
by  iterating  and  converging  on  a  Final  solution  in  the  sense  of  the  example  just  quoted.  The  iterative 
approach  performed  in  this  thesis  is  valid  for  steady-state  problems,  however  the  general  boundary  element 
technique  has  the  capability  to  model  non-steady  problems.  The  problem  might  just  as  easily  been 
formulated  as  the  long-term  response  of  a  non-steady  fluid/structure  system. 

The  major  thrust  of  this  thesis  is  to  research  the  application  of  the  boundary  element  technique  to 


compute  the  incremental  deformation  of  an  elastic  structure  in  fluid  flow.  The  boundary  element,  or 


boundary  integral,  method  has  been  used  to  analyze  the  structure  -  fluid  interaction  problems  in  water 
waves.  As  described  in  References  6  and  11a  semi-infinite  (half)  domain  was  formed  using  the  surface  of 
the  ocean  as  a  free  boundary,  and  the  ocean  bottom  as  a  solid  boundary.  The  solution  was  simultaneously 
solved  using  an  unsteady  approach  with  frequency  domain  decomposition. 

The  boundary  element  technique  was  chosen  over  finite  differences  and  finite  elements  because  of  its 
ability  to  numerically  approximate  both  the  fluid  and  elastic  structure  governing  equations  with  a  common 
definition  of  the  fluid/structure  boundary. 

In  the  iterative  process  the  flow  conditions  (velocity  and  density  of  the  fluid  at  infinity),  the 
structural  conditions  (Young's  modulus  and  Poisson's  ratio),  and  the  original  boundary  geometry  specify 
the  critical  state.  The  first  step  of  the  process  is  to  compute  the  fluid  pressure  field  about  the  structure,  and 
then  determine  the  structural  deformations  caused  by  the  pressure  distribution.  The  new  boundary  geometry 
is  automatically  calculated  from  the  deformations,  and  the  pressure  field  recalculated.  This  process  can  be 
iterated  to  show  how  the  structure  changes  shape. 

An  application  of  this  program  technique  can  be  seen  in  a  flight  test  program  currently  being 
performed  at  the  4950th  Test  Wing,  Wright-Patterson  Air  Force  Base.  The  MILSTAR  airworthiness 
program  is  a  flight  test  program  to  certify  the  airworthiness  of  a  ra dome  to  cover  a  satellite  antenna  for  the 
MILSTAR  terminal  system  (4).  Because  of  a  clearance  of  only  one-half  inch  between  the  antenna  and  the 
radome,  one  objective  of  the  flight  test  program  was  to  measure  radome  deformations  to  verify  adequate 
clearance.  The  analytical  prediction  of  the  radome  deformations  became  a  major  task  of  the  test  team,  and 
was  found  to  be  highly  manpower  intensive.  The  aircraft/radome  geometry  was  first  entered  into  a 
computational  fluid  dynamics  program  called  QUADPAN  (10)  which  calculated  the  pressures  about  the 
radome  assuming  suW  •  steady-state  flow.  QUADPAN  itself  is  a  type  of  boundary  element  code.  The 
radome  structure  w^s  also  r  Meled  with  a  NASTRAN  finite  element  geometry  grid.  Both  of  these  tasks 
took  several  week,  o  set  u  The  users  then  ran  the  QUADPAN  code  at  the  desired  flight  conditions,  fed 
the  resulting  pressures  through  a  program  to  translate  them  from  the  center  of  the  QUADPAN  panels  to  the 
NASTRAN  nodes,  then  input  the  translated  pressure  into  the  NASTRAN  code  to  result  in  the  final  radome 
deformation.  Information  on  the  changes  to  the  pressure  field  due  to  the  deformations,  and  the  resulting 
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change  in  deformation  would  require  extensive  remodeling  to  the  QUADPAN  and  NASTRAN  geometry 
models  after  each  iteration.  Due  to  the  extensive  effort  involved  in  this  remodeling,  no  iterations  were 
performed  on  the  pressure/deformation  solutions.  Deformations  recorded  during  flight  were  used  to 
determine  the  accuracy  of  the  predicted  deformations. 

The  boundary  element  code  presented  in  this  thesis  has  several  advantages  to  the  above  method.  The 
users  would  only  be  required  to  set  up  the  original  aircraft  and  radome  boundaries  which  is  quite  similar  in 
work  effort  to  the  generation  of  the  QUADPAN  model.  Once  the  flow  conditions  and  structural  properties 
were  input,  the  program  would  then  calculate  the  radome  deformations  and  iterate  to  a  steady-state  solution. 
A  more  extensive  unsteady  boundary  element  program  than  that  presented  in  this  research  could  conceivably 
predict  radome  flutter  using  a  time-marching  scheme. 

Subsequent  chapters  of  this  thesis  covers  a  brief  theoretical  background  of  the  boundary  element 
models  used  in  the  computer  program,  a  description  of  the  implementation  of  the  theory,  results  of  the 
investigation,  and  conclusions  and  recommendations  drawn  from  the  research  performed.  While 
simplifying  assumptions  have  been  made  in  the  execution  of  this  investigation,  the  general  boundary 
element  technique  has  the  ability  to  model  complex  governing  equations  and  higher  order  approximation 
techniques.  The  initial  results  obtained  during  the  course  of  this  work  show  promise  towards  the 
suitability  of  the  boundary  element  method  for  such  as  task,  and  follow-on  studies  are  recommended  to 
further  develop  a  code  usable  in  the  field. 
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(2:  Chap  2,  3:  Chap  1) 


During  the  recent  decades  the  emphasis  of  engineering  computations  has  shifted  from  analytical 
approximating  techniques  to  massive  numerical  approximations.  Some  of  the  most  popular  numerical 
approximation  techniques  today  are  the  finite  difference  method  (FDM),  the  finite  element  method  (FEM) 
and  the  boundary  element  method  (BEM). 

The  first  two  methods  are  'domain'  methods. 

These  techniques  discretize  the  domain  of  the  problem  under  consideration  into  a  number 
of  elements  or  cells.  The  governing  equations  of  the  problem  are  then  approximated  over  the 
region  by  functions  which  fully  or  partially  satisfy  the  boundary  conditions.  (2:1) 

In  the  finite  difference  method  a  series  of  nodes  in  a  domain  grid  is  defined  and  the  discretized  version  of  the 

governing  equation  is  satisfied  only  at  those  nodes.  In  the  finite  element  method,  the  governing  equation  is 

satisfied  in  an  average  sense  over  a  region  or  element.  In  the  case  of  the  FEM  the  integrations  are 

conducted  over  the  domain.  In  both  methods,  the  user  must  discretize  the  domain  as  well  as  the  boundary 

of  the  region  under  consideration. 

The  third  method  is  known  as  the  boundary  integral  method,  or  the  boundary  element  method,  and 
satisfies  the  governing  equations  throughout  the  domain  but  only  approximates  the  boundary  conditions. 
One  of  the  first  engineering  applications  of  the  BEM  was  the  solution  for  Laplace-type  problems  by  O.D. 
Kellogg  in  1953  (7).  J.L.  Hess  was  also  an  important  contributor  to  the  development  and  use  of  boundary 
element  methods  in  fluid  mechanics  problems  (5).  Integral  equation  techniques  have  been  in  development 
since  the  early  1950s,  but  they  were  overshadowed  by  the  more  popular  finite  difference  and  finite  element 
techniques  due  to  the  difficulty  of  defining  the  appropriate  Green's  functions  for  the  BEM.  Recently, 
engineers  have  begun  to  rediscover  boundary  element  methods  and  their  advantages.  A  few  of  these 
advantages  are; 

a.  The  governing  equations  are  reduced  from  domain  integrals  to  boundary  integrals  thereby  reducing 
the  order  of  the  problem  by  one  spatial  dimension. 
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b.  Boundary  method  techniques  can  be  coupled  with  other  techniques  to  improve  the  accuracy  of  the 


solution. 


c.  The  boundary  method  can  simply  model  problems  with  infinite  domains. 

<L  Finally,  and  most  important  to  the  application  of  this  thesis,  only  the  boundary  of  the  domain 
must  be  discretized.  Once  the  boundary  geometry  has  been  defined,  both  interior  and  exterior  problems  may 
be  solved  without  discretizing  either  domain.  This  advantage  facilitates  the  modeling  and  tracking  of 
moving  boundaries. 

For  practical  engineering  problems  (full  aircraft,  three-dimensional  geometry)  current  computational 
ability  favors  a  mixture  of  boundary  integral  formulations  where  applicable,  augmented  by  finite  difference 
methods  for  regions  where  non-linear  effects  dominate. 

All  three  of  the  approximation  techniques  are  numerically  related,  and  may  be  derived  through  the  use 
of  the  weighted  residual  technique.  As  an  example  of  the  weighted  residual  formulation  consider  a 
potential  function  u  on  a  domain  Q  which  satisfies  the  governing  equation: 


V  u  -  b  =  0  in  domain  Q 
The  boundary  conditions  are  of  two  types 


(1) 


(a)  Essential  or  Dirichlet  conditions,  such  as  u  =  u  on  f  j 


(b)  Natural  or  Neumann  conditions,  such  as 


du 

dn 


du  _  r 
_ .  cn  r2 


a) 

(3) 


The  barred  terms  are  the  known  conditions.  The  entire  boundary  is  composed  of  Tj  and  r2.  The  outward 
normal  is  defined  as  n  (Figure  1). 

When  the  function  u  and  its  boundary  conditions  are  numerically  approximated,  errors  are  introduced 
into  the  problem  solution.  A  weighting  function  w*  with  continuous  first  derivatives  is  used  to  distribute 
the  error  over  the  domain  and  boundary. 
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The  weighted  residual  statement  becomes; 


|  (  V  u  ~  b)  w  dil  =  0 


in  domain  Q 


(4) 


Further  discussion  on  the  weighted  residual  technique  is  found  in  Appendix  C.  The  finite  difference  method 
can  be  obtained  from  Equation  4,  and  is  called  the  "original  statement"; 


For  finite  differences,  the  basis  functions  u  and  w*  are  normally  different,  where  w*  is  usually  the  Dirac 
delta  function  and  the  derivatives  of  u  are  local  expansions  such  as  Taylor  series. 
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Integrating  Equation  4  oy  pans  once  results  in  the  "weak  statement"; 


Finite  element  techniques  base  numerical  approximations  on  this  form  of  the  governing  equation.  Most 
finite  element  models  are  based  on  the  method  of  Galerkin  where  both  u  and  w*  are  related  so  that  w*  is  a 
variation  of  u  associated  with  virtual  displacements  or  velocities. 

Integrating  Equation  6  by  parts  once  more  obtains  the  "inverse  statement”; 


(7) 


Equation  7  can  be  rewritten  as; 


Equation  8  is  the  basis  for  the  boundary  element  method.  Generally  the  weighting  function  is  chosen  so 
that  the  left-side  of  the  equation  vanishes  leaving  the  boundary  terms. 


Boundary  Element  Formulation  (2:  Chap  3, 3:  Section  2.4) 
If  the  boundary  conditions  are  now  defined  as: 


du 

q  =  ^  and  q 


dw_* 

an 


(9) 
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the  boundary  element  formulation  is: 


|  (v  w  )ud £2=-  f  qw  dT -  f  q  w 

•V  Jr2  Jr , 


dT  + 


I  uqdT  +  I  uqdT  +  I  bw 

yr2  'r,  ■'fl 


dQ  (10) 


Assume  that  a  concentrated  charge  is  acting  at  point  'i'.  The  governing  equation  becomes: 

2  i 

V  u  +  A  =  0  (ID 

where  A‘  is  a  Dirac  delta  function  at  point  T.  The  solution  to  this  equation  is  called  the  fundamental 
solution  where  the  effects  of  boundaries  at  infinity  are  considered.  If  equation  11  is  satisfied  with  the 
fundamental  solution  then: 


f  2  . 

I  u(  V  W 

JQ. 


)d£2  =  u 


(12) 


Equation  10  now  becomes: 


u+|  uqdT+l  u  q  dT  =  I  qudT+l  qudT+l  bw 

Jr2  Jrl  Jr2  Jr1  Jn 


dn 


(13) 


Equation  13  is  valid  for  any  point  in  domain  £2,  but  must  be  taken  to  the  boundary  to  solve  for  the 
boundary  conditions.  To  do  this  a  region  with  radius  e  is  formed  around  the  boundary  point  i  (Figure  2). 
Equation  13  is  separated  into  integrals  around  Te  and  around  r-re.  By  taking  the  limit  as  e  approaches  zero, 

the  integrals  over  T  -  Te  are  continuous.  The  integral  over  Te  for  the  right  hand  side  contains  a  weak 
singularity  but  remains  continuous.  The  left  side  integral  contains  a  higher  order  singularity  because  of 
the  derivative  term.  As  the  point  i  approaches  the  boundary  along  the  normal,  it  produces  a  discontinuity 
as  the  point  passes  through  the  boundary. 
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Figure  2.  Applying  the  BEM  Equation  to  the  Boundary 


On  the  boundary  Equation  13  is  reduced  to: 


J  uq  dT  +  J 

f  u  q*  dT  =  J 

\  qu*dT  +  J 

f  qu  dT  +  J  bw  dQ 

(14) 

Jr2  J 

r2 

The  term  (2  rc  c‘)  represents  the  discontinuity  and  is  a  function  of  the  boundary  shape.  By  employing  a 
method  similar  to  that  for  evaluating  the  discontinuity,  it  can  be  shown  that  the  term  equals  the  internal 
angle  of  the  boundary  at  the  point  i  (3:  63,64).  For  a  smooth  boundary  the  term  equals  n  and  c1  becomes 
1/2.  A  second  method  of  calculating  c*  will  be  discussed  in  Chapter  3. 

Direct  and  Indirect  Methods 

The  direct  formulation  of  the  BEM  is  the  formulation  discussed  above.  It  may  be  derived  through  the 
use  of  we.0hted  residuals  as  shown,  or  through  Green's  third  identity,  Betti's  or  similar  theorys.  The  final 
solution  is  the  function  u  and  its  derivatives.  For  example,  when  modeling  the  elastic  structure  the  final 
solution  is  the  deformation  at  the  node  points.  An  advantage  of  the  direct  formulation  is  that  a  general 
surface  with  comers  or  edges  may  be  used. 

The  indirect  formulation  of  the  Laplace  equation  problem  sets  up  the  problem  as  solely  a  single-layer 
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potential  or  solely  a  double-layer  potential  of  continuous  source  distributions  over  the  boundary.  The 
indirect  method  carries  a  restriction  that  the  boundary  must  be  smooth.  The  equations  are  then  solved  for 
the  source  densities.  Further  discussion  on  the  indirect  method  may  be  found  in  Reference  7. 

Elements 

As  the  name  'boundary  element  method’  implies,  the  boundary  of  domain  Q  is  divided  into  segments 
called  elements  (Figure  3).  The  general  theory  places  no  restrictions  on  the  modeling  of  the  element.  For 
simplicity,  linear  (flat)  elements  were  used  in  the  development  of  this  thesis,  but  higher  order  elements 
may  be  used  to  model  the  boundary  more  smoothly. 


Element 


b)  4  quadratic  elements  (8  nodes) 


Figure  3.  Geometry  Elements 

The  order  of  the  geometry  model  will  determine  the  continuity  of  the  geometry  slope.  Throughout  a 
continuous  boundary  the  du/dn  terms  may  be  calculated  everywhere,  but  on  a  discontinuous  boundary  the 
du/dn  terms  are  undefined  at  the  ends  of  the  elements  because  the  normal  is  undefined.  A  discontinuous 
geometry  slope  will  also  cause  the  solution  to  display  sharply  singular  behavior  in  the  vicinity  of  the 
comers.  The  major  disadvantage  to  continuous  geometry  models  is  the  degree  of  effort  involved  in  the 
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implementation.  For  the  purposes  of  this  investigation,  linear  geometry  modeling  was  used  and  therefore 
the  geometry  slope  is  discontinuous  at  the  element  comers.  The  consequences  of  using  discontinuous 
geometry  are  discussed  in  Chapters  4  and  5 

Approximations  of  the  function  u  and  its  derivatives  are  assumed  to  vary  over  the  element.  The 
governing  equation  is  then  written  at  various  locations  on  each  element  called  nodes.  The  number  of  nodes 
on  each  element  corresponds  to  the  order  of  the  approximating  functions  of  u  and  q,  and  not  necessarily  to 
the  order  of  the  element  geometry.  For  example,  u  may  be  modeled  with  a  parabolic  function  ,  and 
therefore  three  function  nodes,  while  the  geometry  element  may  be  linear  (Figure  4).  The  functions  u  and  q 
do  not  have  to  be  modeled  with  the  same  order,  in  fact  it  may  be  more  advantageous  to  take  q  of  one  order 
less  than  u  because  it  is  the  derivative  of  u.  For  the  purpose  of  this  research,  the  approximations  of  u  and  q 
were  equivalent  order  for  simplicity  of  implementation. 

G - U - Q - V - Q - U - G 


Order  of  Geometry  =  1st  (2  nodes) 

Order  of  Function  u  =  2nd  (3  nodes) 

Order  of  Function  q  =  1st  (2  nodes) 

Figure  4.  Approximating  Order  of  Functions  and  Corresponding  Nodes 

Internal  and  External  Domains 

At  each  node  an  outward  normal  to  the  element  is  calculated.  The  outward  direction  is  defined  in 
relation  to  whether  an  internal  or  external  problem  is  being  solved.  An  internal  problem  is  one  which  the 
domain  under  study,  ft,  is  enclosed  by  the  boundary,  I~.  When  the  domain  is  outside  the  boundary,  it  is  an 
external  problem  (Figure  5).  In  the  case  of  potential  flow  past  an  elastic  structure,  the  flow  would  be  an 
external  problem  and  the  structure  is  an  internal  problem.  The  ability  of  the  boundary  element  technique  to 
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solve  both  internal  and  external  solutions  with  the  same  boundary  geometry  allows  the  user  to  easily  define 
and  track  moving  boundaries. 


a)  Internal  Problem  b)  External  Problem 

Figure  5.  Internal  and  External  Problem  Definitions 
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The  solution  technique  first  calculates  the  potential  of  the  fluid  about  the  structure.  Next,  the 
velocities  and  pressures  on  the  surface  due  to  the  potential  are  determined.  Finally,  the  elastic  deformations 
from  the  pressure  field  are  found.  A  direct  boundary  element  formulation  was  used  for  the  potential, 
velocity  and  deformation  models.  This  research  was  limited  to  a  two-dimensional  problem,  and  the 
following  model  development  will  only  discuss  two-dimensional  solutions.  A  three-dimensional  problem 
is  conceptually  solved  in  a  similar  manna*. 

Potential  Flow 

The  beginning  equation  of  the  direct  method  boundary  element  formulation  for  potential  flow  is: 

Ci©Ui(|)  +  J  u(x)  q  (^)  dT(x)  a  J\j(x)  w  fox)  dT(x)  (15) 

where; 

£  are  the  locations  of  a  continuous  distribution  of  sources  throughout  domain  ft 
x  are  the  coordinates  for  the  jtoint  where  the  potential  is  being  evaluated 
T  is  the  boundary  of  the  domain 

i  is  a  point  where  the  fundamental  solution  is  applied,  and  is  called  a  nodal  point 
u(x)  is  the  unknown  and  is  the  potential  at  point  x 
q  =  du/dn 

w*(£,x)  is  the  weighted  residual  weighting  function 
q*(^,x)  =  dw*/dn 

The  fundamental  solution  to  Laplace's  equations  is  used  as  the  weighting  function  w*.  This  solution 


is  the  free  space  Green's  function  for  the  governing  equation. 


For  two-dimensional  problems  the  solution  is: 


Equation  15  has  been  divided  through  by  2  n  so  that  the  2  n  appears  in  the  denominator  of  w*  and  q' 
distance  between  £  and  x  is  shown  in  Figure  6  and  is  determined  by: 


=  15 -x|  =  V(x  <£)  -  x  (x))2  + 


(y  ft)  -  y  W)2 


for  a  Cartesian  two-dimensional  system. 


Figure  6.  Source  and  Field  Points 


(16) 


.  The 


(17) 


14 


When  the  boundary  is  discretized  into  elements  Equation  15  becomes: 


N  C  *  N  C  * 

ci  Uj  +  £  J  u  q  dT  =  £  J  w  qdT  (18) 

j  =  l  ri  j  =  l  ri 

where  N  is  the  number  of  geometric  elements  and  Tj  is  the  length  of  element  j. 

An  elemental  coordinate  system  is  set  up  for  each  element  as  shown  in  Figure  7.  The  values  of  u 
and  q  vary  over  the  element  as  a  function  of  the  coordinate  tj. 


Figure  7.  Elemental  Coordinate  System 


The  functions  u  and  q  may  be  approximated  with  shape,  or  interpolation,  functions  such  that: 

i U| ) 

U  0D=  N1(ti)u1+  N2(o>2+  .  .  .  +  Nr(n)ur  =  [N1(Ti),N2('n),  .  .  .  Nr(n)  ]  U2  j  (19) 
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q(n>=  Nl(TOq,  +  N2(Ti)q2+  .  .  .  +  Nr(n)qr  =  [n,(tO,N2(ti).  .  .  .  Nr(n)  ] 


qi 
< b 

qr 


(20) 


where  r  is  the  order  of  the  function  modeling  (i.e.,  second,  third,etc).  Appendix  C  provides  a  brief  overview 
of  shape  functions.  Notice  that  ur  and  qj.  are  the  amplitudes  of  the  functions  at  the  node  points  and  are  no 
longer  functions  of  T|.  Therefore,  they  may  be  taken  outside  the  integrals. 

The  integrals  are  with  respect  to  the  boundary  r  ,  but  u  and  q  are  functions  of  q.  A  Jacobian  is  used 
to  transform  the  integral  as  detailed  in  Appendix  C  ; 

cG  =  |G|dn  (21) 

By  using  Equations  19, 20  and  21,  Equation  18  now  becomes: 


Over  each  element  j  there  are  r  nodes.  Globally,  there  are  m  =  j  x  r  total  nodes  on  the  boundary  (Fig  8). 
The  integral  on  the  left-hand  side  of  Equation  22  can  be  written  as: 


A 
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where  the  h  jj  are  influence  coefficients  defining  the  interaction  between  the  point  i  under 

consideration  and  a  particular  node  k  on  element  j.  The  same  can  be  done  for  the  right-hand  side  integral, 
so  that  Equation  22  can  now  be  written  as: 


N  'T  N  k 

CiUi  +  I  hijUk  »  J  g  ijQk 

j  =  l  j  =  l 


(24) 


Figure  8.  Elemental  and  Global  Node  Numbering 
When  summed  over  all  the  elements,  equation  24  can  be  written; 


where  m  equals  the  total  number  of  nodes  around  the  boundary  in  the  global  numbering  system 
and: 


Hi, 


>/n„, 

J  r 


j(n)q(t*)lG|dn 


(26) 
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Gim  =  /  Nnj  (TO  w  &x)  |G|  dq  (27) 

*  r 

where  the  index  i  is  on  the  node  i  and  nj  is  the  nth  node  on  element  j  in  the  global  numbering  system. 
Using  Figure  8  for  reference,  Hj  15  is  the  influence  at  global  node  15  ( n  =  3,  j  =  5)  due 


to  point  1.  The  qq  term  may  be  aritted  to  the  H  matrix  such  that 

Hij  =JHij  fori*j 
Hjj  =  Hjj+  Cj  for  i=j 

Equation  25  is  repeated  for  each  node  i  (m  total)  and  produces  a  matrix  equation; 

H  u=  G  q  (29) 

For  a  smooth  boundary,  it  can  be  shown  that  the  coefficient,  q,  is  equal  to  1/2  as  discussed  in  Chapter  2. 
If  the  surface  is  not  smooth,  it  will  not  be  equal  to  1/2.  The  diagonal  term  of  H  may  also  be  calculated  by 
using  the  fact  that  when  a  uniform  potential  is  applied  on  the  whole  body  the  sum  of  the  values  of  q’s  must 
be  zero.  For  internal  or  closed  problems ; 

HI  =  0  (30) 

so  that  the  sum  of  each  row  in  H  should  be  zero  (3:  111).  Therefore, 


Hi  j  =  -  £  Hjj  i  *  1,2 . 

j  =  l 

For  exterior,  or  unbounded,  problems  it  can  be  shown  (3: 1 12)  that: 


(31) 


Hu  =  -  I  Hij  +  1  i  =  1,  2, . N,  i  *  j  (32) 

j  =  l 

The  boundary  conditions,  q,  can  be  calculated  using  the  assumption  of  zero  flow  through  the  boundary. 
This  means  that  the  source  must  produce  equal  and  opposite  flow  to  the  flow  at  the  boundary,  or: 

=-n>'Vom  (33) 
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where  nj  is  the  outward  unit  normal  at  point  i  and  VM  is  the  flow  condition  at  infinity. 

The  matrices  H  and  G  may  be  calculated  numerically  using  numerical  quadrature  for  the  integral 
terms.  Equation  29  is  a  m  x  m  system  of  equations.  As  the  integration  progresses  around  the  boundary, 
for  each  node  i  there  will  be  a  panel  for  which  the  integration  is  singular  (r  approaches  and  equals  0).  Care 
must  be  taken  to  properly  integrate  these  singular  panels.  Note  that ; 


q'<«  = 


9n  3n  \2«  r 


2n 


t  *  V 1  In  — 
r 


On  a  flat  singular  panel  the  normal  is  perpendicular  to  the  panel  and  therefore; 


=  0 


(34) 


(35) 


Therefore,  the  H  terms  are  zero  on  the  singular  panel.  With  higher  order  geometric  modeling  this  is  not 
the  case  and  this  term  must  be  addressed. 


Velocity  and  Pressure  Distribution 

Once  Equation  29  is  solved  for  the  potential  at  the  nodes,  the  velocity  and  pressure  distribution  at  the 
nodes  can  be  determined.  The  components  of  the  perturbation  velocity  at  node  i  are  calculated  using: 


A  m  x  m  system  of  equations  may  be  set  up  for  Equation  36  similarly  to  that  described  above,  and 
numerical  quadrature  used  to  calculate  the  terms  of  the  matrices.  The  functions  q(x)  and  u(x)  are  known 
from  previous  calculations  so  that  the  velocity  may  be  solved.  The  components  of  the  perturbation 
velocity  in  the  x  and  y  directions  are: 


u  -|j&  and  v  -^§1 
ox(x)  dy(x) 


(37) 
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and  the  total  perturbation  velocity  is  calculated  by: 


Vp  =  Vu+v  (38) 

The  total  velocity  at  node  i  is  equal  to: 

V  =  V..-Vp  (39) 

The  pressure  field  about  the  body  is  computed  using: 


P  =  ip.V-Cp  (41) 

Px=P'nx  and  Py=P'ny  (42) 

Px  and  Py  are  the  components  of  the  pressure  acting  on  the  surface  or  the  surface  tractions. 


Elastic .  S  tructure 


The  development  of  the  boundary  element  formulation  for  the  elastic  structure  is  quite  similar  to  that 
of  the  potential  flow.  The  structure  consists  of  subregions  formed  by  an  outer  boundary,  inner  boundary 
and  connecting  elements  as  shown  in  Figure  9.  The  connecting  elements  allow  the  effect  of  a  force  at  one 
point  on  the  boundary  to  be  transferred  around  the  structure.  For  this  investigation,  the  same  number  of 
nodes  were  used  on  the  outside,  inside  and  connecting  elements  although  this  is  not  required. 

The  beginning  statement  for  each  subregion  is: 


CijUj  + 


J  PijUjdT  =  J  UijPjdT  +  J  UijbjdQ 

*  r  'r  •  n 


(43) 


where 
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uj  is  now  the  unknown  displacements  in  the  j  direction 

Pj  is  the  traction  or  natural  boundary  conditions  in  the  j  direction  (the  pressures  calculated  above) 
bj  is  the  body  force  in  the  j  direction  and  is  assumed  zero  for  this  research. 


The  fundamental  solution  w*  for  a  two-dimensional  plane  strain  problem  is  (2: 187): 


w*j  = - l- -  1(3  -  4v)  ln(r)  8{i +£-£-) 

8  re  (1-  v)  G  l  axidXjj 


(1-  v) 

The  solution  p*  is  defined  as: 


(44) 


* 

Pij  = 


-1 


4  re  (1-v)  r 


(l-2v)5ij  +2|L|L 

dXj  oXjJ 


|l-(I-2v)(^nr|Lni) 
on  ox  j  ox: 


(45) 
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As  with  the  potential  flow  equations,  the  functions  u  and  q  are  represented  with  shape  functions  and 
the  integrals  transformed  from  the  dT  to  the  dq  system. 

Equation  43  can  be  written: 


m  is  the  total  number  of  nodes  around  the  outer  boundary  and  is  equal  to  the  nodes  around  the  inner 
boundary  and  also  around  the  connecting  elements.  The  approximating  functions  of  u  and  q  are  modeled 
with  the  same  order,  r.  on  all  elements.  Therefore,  m  =  3  x  r  nodes  around  the  region.  Note  that  the 
matrix  terms  are  2  x  2  for  the  two-dimensional  case  because  they  contain  the  terms; 


^ij  xx  hy  ^ 
**ijyx  kijyjr 


8ij  = 


8  ij  m  8  ij 
8yj*  8ijjy 


(47) 


For  the  three-dimensional  case  the  matrices  would  be  3  x  3. 

Equation  46  produces  a  lx(3xrx2)  system  of  equations  for  each  node  i.  Equation  46  is  calculated  twice 
for  each  of  the  nodes  on  the  connecting  element;  once  for  the  element  to  the  left  and  once  for  the  element 
to  the  right.  In  this  way  the  matrices  have  (4xrx2)  rows.  The  final  system  of  equastions  is  a  (4xrx2)  x 
(3xrx2)  system. 

The  coefficient  Cj  may  be  handled  similarly  to  Equation  29  to  form  the  system  of  equations: 

Hu  r  GP  (48) 

The  Cj  coefficient  may  be  calculated  using  rigid-body  considerations  in  a  bounded  body  (3:201).  If  a  unit 
rigid-body  displacement  is  made  in  any  direction  Equation  48  becomes: 

H  I  d  =  0  (49) 
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where  I  ^  is  a  vector  of  unit  displacements  in  the  d  direction.  The  diagonal  terms  of  H  become: 


N 

l*i  i  ~  "  2  h  i  j 
1  =  1 


fori  *  j 


(50) 


Assembly  of  Subregions.  For  each  subregion  in  the  structure  the  system  in  Equation  48  may  be 
developed.  To  solve  for  the  deformations  about  the  entire  body,  the  systems  of  equations  must  be 
assembled  into  a  global  system  of  equations.  This  is  done  by  renumbering  the  nodes  around  the  structure. 
The  nodes  around  the  outside  are  numbered  1  to  m,  the  nodes  about  the  inside  m+1  to  2m,  and  the  nodes  on 
the  connecting  elements  2m+l  to  3m  (Figure  10).  The  H  and  G  matrices  are  now  (8xm)  x  (6xm).  The 
global  system  of  equations  is: 
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where  H'jik  and  G^  are  the  influence  coefficients  at  node  i  due  to  panel  j  in  the  Ik  direction.The  boundary 
forces  around  the  boundaries  are  the  pressures  from  the  fluid  and  have  been  calculated  previously.  The 
boundary  tractions  along  the  connecting  surfaces  are  unknowns  and  must  be  calculated  with  the  unknown 
deformations.  The  columns  of  G  and  the  corresponding  rows  of  P  are  transferred  to  the  left-hand  side  of 
the  equations  forming  a  (8xm)  x  (8xm)  matrix  with  (8xm)  unknowns.  This  final  set  of  equations  may  now 
be  solved  for  the  unknown  deformations  and  internal  tractions. 


IV.  Implementation 


Once  the  potential,  velocity  and  elastic  deformation  models  were  developed,  they  were  coded  into  a 
Fortran  computer  program.  Several  simplifications,  assumptions  and  numerical  approximations  were 
included  in  the  program  code. 


Geometry 


The  structure  used  during  this  research  was  a  two-dimensional  thick-walled  cylinder  with  an  outer  and 
inner  radius  (Figure  1 1).  The  geometry  was  modeled  as  flat  elements  with  geometry  nodes  at  the  left  and 
right  ends  of  the  elements  (n=-l,  tj=1).  The  number  of  elements  around  the  cylinder  could  be  varied  so 
that  the  cylinder  could  be  modeled  very  crudely  with  8  elements  as  shown  in  Figure  11,  or  finely  with  up 
to  150  elements.  The  number  of  outer  boundary  elements,  inner  boundary  elements  and  connecting 
elements  were  equal  for  ease  of  geometry  definition.  The  elements  were  numbered  in  a  counterclockwise 
direction  with  the  structure  domain  remaining  on  the  left-side  as  one  travels  around  the  outer  boundary,  and 
on  the  right-side  along  the  inner  boundary.  The  connecting  elements  were  specified  so  that  the  left  end 
would  be  on  the  outer  boundary  and  the  right  end  on  the  inner  boundary. 

The  element  Jacobian  was  calculated  as  detailed  in  Appendix  C,  and  the  outward  normals  calculated 
by: 


Jy/ 


n,  =- 


'*1 


Jacobian 


and  nv 


Jacobian 


(52) 


Because  of  the  counterclockwise  setup  of  the  elements,  the  outward  normal  on  the  outer  boundary  is 
into  the  structure  and  therefore  the  correct  sign  for  the  potential  flow  calculations.  For  the  elastic 
deformation  calculations  a  negative  sign  must  be  applied  to  the  normal  so  that  it  may  be  an  outward 
normal  for  the  structure.  The  inner  boundary  normal  is  in  the  correct  direction  for  the  elastic  deformations. 
The  correct  direction  of  the  connecting  elements  outward  normals  were  determined  in  the  same  manner. 


26 


Figure  1 1.  Geometry  Modeled  With  8  Flat  Elements 


The  order  of  the  potential,  pressure  and  deformation  approximation  was  equal  and  could  be  specified 
from  constant  (1  node)  to  4th  order  (5  nodes).  Therefore,  the  potential,  pressure  and  deformation  were  all 
calculated  at  the  same  node.  The  general  BEM  formulation  does  not  require  this  assumption.  The  distance 
of  the  end  nodes  to  the  element  comer,  S,  was  specified  by  the  user.  In  this  way  the  effect  of  the  distance 
between  the  end  node  and  the  comer  could  be  investigated.  The  end  node  was  never  placed  directly  on  the 
comer  (8  =  0).  As  previously  discussed  there  are  several  difficulties  that  arise  if  the  node  is  placed  directly 
on  the  comer  which  is  a  geometry  slope  discontinuity.  The  remaining  nodes  were  spaced  equally  between 
the  end  nodes  (Figure  12)  with  the  same  number  of  nodes  used  on  the  outer,  inner  and  connecting  boundary 
elements. 
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Figure  12.  Nodal  Spacing  Within  An  Element 

Once  the  deformations  at  the  node  points  were  calculated,  the  deformations  at  the  ends  of  each  panel 
were  determined  using  the  shape  functions.  For  example,  the  deformations  at  the  left  end  of  a  panel  are: 

ux(ri  =  -l)  =  N,(-l)uxl  +  N2(-l)ux2  +  ...  +  Nr(-l)u„  (53) 

Uy(Ti  =  -1)  =  N,(-l)  uyl  +  N2(-l)uy2+  ...  +  N,(-l)  uyr 

where  Nj(-l)  is  the  shape  function  afq=-l  for  element  node  i,  and  u  ^  and  u  y,  are  the  deformations  at  node 

i  on  the  element.  The  deformations  at  the  right-end  of  the  element  were  calculated  in  a  similar  manner. 

Because  the  discontinuous  slope  at  the  comers  causes  stress  concentrationsat  the  comers  the  deformations 

are  not  continuous  from  element  to  element  In  order  to  force  continuity  of  the  new  geometry,  the 

deformations  as  each  comer  were  averaged  with  the  deformation  of  the  right  panel  and  the  deformation  of 

the  left  panel.  The  averaged  deformations  were  applied  to  the  original  geometry  resulting  in  the  new 

geometry  boundaries.  Once  the  new  outer  and  inner  geometry  was  calculated,  the  connecting  element 

geometry  was  recalculated. 

Potential  How 

Potential  flow  acts  only  on  the  outer  boundary.  There  is  no  flow  inside  the  cylinder.  The  potential 
flow  subroutine  consists  of  an .  outer  loop  on  the  number  of  nodes  and  an  inner  loop  on  the  number  of 
elements.  For  each  node  the  routine  sets  up  Equation  15  and  numerically  integrates  over  all  the  elements. 
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The  functions  w*  and  q*  are  defined  as: 


On  the  non-singular  elements  the  integration  was  performed  with  a  Gaussian  quadrature: 


f 1  f  (ti)  dn  =  I  fttOwi  (55) 

j-i  i=i 

where  N  is  the  order  of  the  quadrature,  T|j  are  the  quadrature  values  and  wj  the  quadrature  weights.  The  order 
of  the  quadrature  used  was  determined  by  comparing  the  numerical  integration  results  with  analytical  results 
as  the  element  end  points  were  moved  towards  the  comer.  In  this  way  it  was  found  that  a  12th  order 
quadrature  was  required  to  obtain  accurate  integrations  for  small  values  of  8.  This  was  also  verified  by 
comparison  of  the  diagonal  terms  of  the  H  matrix  as  the  quadrature  order  was  increased.  Because  the 
geometry  was  always  smooth  at  the  nodes,  the  diagonal  term  should  be  .5.  A  12th  order  quadrature  was 
required  to  obtain  .5  as  the  end  nodes  came  very  close  to  the  comers  (5  =  .01).  This  high  order  is  due  to  the 
1/r  term  which  causes  r  to  become  very  small  as  the  node  is  moved  to  the  end  of  the  panel. 

The  singular  element  integrations  were  performed  with  a  different  quadrature  formula  The  singular 
terms  of  the  H  matrix  (q*  integral)  vanish  because  the  normal  is  perpendicular  to  the  panel  for  flat  panels 
(Eq  34, 35).  For  the  G  matrix,  a  one-dimensional  logarithmic  Gaussian  quadrature  formula  is  used  (3:449). 
Because  the  formula  is  designed  for  an  integral  from  0  to  1  with  the  singularity  at  0,  the  integral  is  broken 
into  two  separate  integrals: 
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r 


-  T1  -  Tl(r  =  0)  = 

■n  =  -  and  ri 

-1  -  Ti(r  =  0) 


1-Ti(r  =  0) 


S  is  now  defined  as  the  location  of  the  singularity  on  the  panel  (ti(r=0)).  On  the  singular  element 


r  =  |S-t||L^L 

where  Length  is  the  length  of  the  panel.  Substituting  Equation  58  into  Equation  56  results  in: 


/  Win}}). 


f  f(T!)ln|lj(l  +  S)dn  -  f  f  (n)  In  f  |l  +  S|^M^j(l  +  S).5dn 
•'O  j -1  ' 


Jf'H; 


+  I  f  (H)  In  l  (l-S)dn  -  f(T1)lnMS-l|i^L  (1  -  S)  .5  dq 


q  =  2q-l  and  q=2q-l  (60) 

Using  Equation  59  the  singular  terms  of  the  G  matrix  were  calculated.  Once  the  H  and  G  matrices  are 
filled,  the  program  calculates  the  diagonal  terms  of  the  H  matrix  using  Equation  32,  and  the  boundary 
condition  matrix  q  using  no  flow  through  the  boundary  (Equation  33). 


The  velocity  and  pressure  subroutine  is  similar  to  the  potential  flow  routine  in  that  it  performs  an  outer 
loop  on  the  number  of  nodes  and  an  inner  loop  on  the  number  of  elements.  The  fundamental  solutions  are 
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and 

dx(x)  2  K  f2  dy(x)  2  7T  r2 


•  2  •  2 

aq  _  1  f  ni-(rini-rini)2ri  8q  _  1  r  n;- (rin1-rini)2  r; 

3x(x)  2  it  4  dy(x)  2  tc  4 


where  the  q  and  rj  terms  are  the  components  of  r  in  the  x  and  y  direction.  The  non-singular  integrals  are 

approximated  by  a  12th  order  Gaussian  quadrature. 

On  the  singular  panels  the  following  approximations  are  made  (Figure  13): 


r  =  x-ti 


Length 


-/x-H 


r  cos  0  and  rf  = 


rsin  0 


Integration  Point  'H 


(x2,y2) 


Singularity  4 


(xl.yl)> 


Figure  13.  Singular  Panel  Approximation  Definitions 
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Because  r  is  always  postive,  a  correction  factor  was  required  to  place  the  correct  sign  on  the  r  components. 
The  angle  formed  by  the  x-axis  and  the  element  is  8  and  is  calculated  by: 


cos  0  =  *2~  —  and 
Length 


sin  0  = 


y2-yi 

Length 


(65) 


where  (*2^2) 318  *e  coordinates  of  the  right  end  of  the  panel,  and  (xj.yj)  are  the  coordinates  for  the  left 
end. 

By  using  Equations  63, 64  and  65  Equation  61  reduces  to  a  1/r  singularity  term.  A  quadrature  using 
a  numerical  evaluation  of  Cauchy  principal  values  (3:451)  was  used  for  the  integrations. 


f  On) 
(n  -s) 


K 

dn  =  Z 

i=l 


f[(-l  -S)tii+  sj  Wj  +  f(S)  In  |-1  -  S  | 


(66) 


*  Sf[d-S)Tli  +  SlWi  -  f  (S)  In  1 1  -  S  | 
H  L  J 


Equation  62  reduces  to  a  1/r^  singularity.  From  Reference  9  this  can  be  approximated  by  a  quadrature  rule 
using; 


I. 


f(n) 


a*,  =  t  .  1^  +  2 


=1 


i=l 


Tli 


f  (T1=0) 


(67) 


Note  that  Equation  67  assumes  that  the  singularity  is  in  the  center  of  the  element.  To  use  this  the 
singular  element  was  divided  into  one  segment  center  around  the  sigularity,  and  a  second  segment  which 
does  not  contain  the  singularity. 


.5  (end  -  start)  dr)  +  non-singular  integral 


(68) 
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where  start  =  start  of  singular  segment  and  end  =  end  of  singular  segment  If  the  singular  segment  is  from 
-1  to  end,  then  the  non-singular  integral  will  be: 


non-singular  integral  = 


,5  (1-  end)  dn 
T1 


where 


2 


(fS) 


- 1 


(69) 


If  the  singular  segment  is  from  start  to  1,  then  the  non-singular  integral  will  be: 


i 


^  »V  #V  /  i 

non-singular  integral  =  .5  (start  +  1)  dr)  where  i\=2  (■  TLt.-L.|  - 1 

~2  \start  +  1  / 


(70) 


The  second  integral  may  be  solved  using  a  normal  Gaussian  quadrature  while  the  singular  segment  may  use 
Equation  67. 

Both  the  boundary  conditions  and  the  potential  are  known  at  the  nodal  points.  Equation  36  is  solved 
at  each  node  in  the  x  direction  for  the  x-component  of  perturbation  velocity  and  in  the  y  direction  for  the 
y-component  of  perturbation  velocity.  The  total  perturbation  velocity  is  subtracted  from  the  oncoming 
flow  velocity  at  each  node  to  determine  the  final  velocity  at  each  node.  The  pressure  field  on  the  cuter 
surface  is  then  calculated  using  Equations  40, 41  and  42.  By  changing  the  density  of  the  fluid,  the  user 
may  increase  or  decrease  the  magnitude  of  the  pressure  on  the  cylinder. 

For  the  pressure  field  on  the  inner  boundary  the  inner  pressure  was  assumed  constant,  and  was  set 
equal  to  a  pressure  on  an  outer  node.  This  approximates  the  cylinder  being  vented  to  the  outside.  The  node 
could  be  chosen  by  the  user,  and  was  selected  as  a  node  on  the  back  of  the  cylinder  where  the  pressure 
approached  stagnation  pressure.  Once  a  vented  node  was  specified,  the  pressure  field  on  the  inner  boundary 
was  determined. 
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The  elastic  deformations  were  determined  in  a  nested  three  step  process.  The  outer  step  is  on  the 
subregions  of  the  structure,  the  next  step  on  the  nodes  around  the  subregion,  and  the  innermost  step  on  the 
elements  around  the  subregion.  The  fundamental  equations,  w*  and  q*,  are  given  in  Equations  44  and  45. 
For  the  non-singular  integrals  the  12th  order  quadrature  is  used.  For  the  singular  integrals  the  equations  are 
reduced  to  ln(l/r)  and  (1/f)  singularities  whose  quadrature  forms  have  been  previously  discussed.  As  the 
routine  progresses  around  a  subregion  it  places  the  terms  of  the  H  and  G  matrices  into  the  global 
numbering  system  rows  and  columns.  When  the  routine  has  finished  looping  around  all  the  subregions, 
the  matrices  are  entirely  filled.  This  routine  did  not  attempt  to  optimize  the  H  and  G  matrices.  If 
optimization  were  to  be  performed,  the  matrices  would  be  banded  which  would  reduce  the  time  required  for 
solution.  The  diagonal  terms  of  the  H  matrix  are  calculated  using  Equation  50  and  the  unknown  tractions 
are  moved  to  the  left-hand  side  of  the  equation.  The  system  of  equations  is  then  solved  for  the  deformations 
at  the  nodal  points. 

Solution  of  Linear  Equations 

The  matrix  system  HU  =  Gq  was  sol  ed  using  a  routine  found  in  Reference  2  (2:67).  The  proper 
execution  of  this  routine  was  verified  using  IMSL  routines.  IMSL  was  not  used  because  it  requires 
additional  workspace.  Because  of  the  large  systems  of  matrix  equations  used  in  the  program,  memory  space 
was  a  driving  factor  to  the  size  of  the  problem  solved. 


Continuity  Considerations 

Because  the  end  nodes  of  the  potential  and  deformation  functions  were  not  on  the  end  of  the  element, 
the  approximating  functions  were  not  constrained  to  be  continuous  from  element  to  element  To  force  the 
potential  to  be  continuous  between  elements  a  continuity  restriction  could  be  enforced  at  the  comers.  At 
each  comer  the  function  u  computed  using  shape  functions  from  the  left  element  equalled  the  function 


computed  using  shape  functions  from  the  right  element: 

Ni(-l)ullt+N2(-l)u2tt+  . . .  +Nr(-l)urrt  =N1(l)unf+N2(l)u21f+  . . .  +Nr(l)urlf  (71) 

on 

Ni(-l)ultt+N2(-l)u2n+  . . .  +NI(-l)uIfl  -NjOJuu-NjOJum-  . . .  -Nr(l)urlf=0  (72) 

A  row  of  the  Hu  =  Gq  system  was  replaced  on  each  element  by  Equation  72.  When  solved,  the  new 
system  of  equations  produced  continuous  values  of  the  function  about  the  surface.  In  practice,  the  order  of 
the  approximating  function  was  increased  to  allow  the  same  number  of  collocation  points  on  the  elements 
as  were  used  without  the  continuity  restraint  An  alternative  approach  not  used  here  would  have  been  to 
use  a  least  squares  solution  by  applying  a  single  value  decomposition  solver.  This  approach  would  not 
have  ensured  continuity,  however,  and  was  rejected  in  favor  of  the  first  approach. 

Slope  continuity  would  be  desirable.  However,  since  the  geometry  was  modeled  as  straight  line 
segments  the  solution  was  expected  to  display  singular  behavior  in  the  vicinity  of  the  comers  (8).  This  is 
analogous  to  the  stress  concentrations  expected  at  sharp  comers  in  elastic  solutions.  If  the  geometry  model 
is  forced  to  have  continuous  slope  in  the  future,  it  would  be  desirable  to  enforce  not  only  continuity  of  the 
unknown  but  also  continuity  of  its  slope  at  the  junction  of  the  elements. 

The  notion  of  a  least  squares  approach  could  be  applied  to  higher  order  elements.  This  could  replace 
the  existing  linear  terms  satisfying  the  original  geometric  nodal  data  and  the  slope  continuity  requirement. 
The  new  geometric  nodal  points  could  then  be  used  to  describe  the  boundary. 

Iterative  Process 

The  potential,  velocity  and  pressure  field  due  to  the  original  geometry  are  calculated  the  first  pass 
through  the  routine.  The  deformations  caused  by  the  pressure  field  (surface  tractions)  are  determined,  and 
the  average  comer  deformations  are  applied  to  the  original  geometry  resulting  in  a  deformed  geometry. 

A  new  potential,  velocity  and  pressure  field  are  calculated.  Beginning  with  the  second  iteration  the 
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elastic  structure  is  in  a  pre-stressed  condition.  If  the  total  pressure  field  is  used  to  calculate  the 
deformations,  it  would  be  equivalent  to  releasing  the  stresses  in  the  structure  after  each  iteration  and 
beginning  from  an  initial  condition.  To  model  the  pre-stressed  condition  the  new  surface  tractions  are 
subtracted  from  the  old  tractions.  This  difference  determines  the  tractions  caused  only  by  the  deformations 
from  the  previous  iteration.  The  difference  is  then  used  to  calculate  the  new  deformations  and  new 
geometry. 

The  program  iterates  in  the  above  manner,  calculating  new  tractions  and  deformations  each  iteration. 
If  the  structure  has  adequate  stiffness,  the  incremental  deformations  will  approach  zero,  and  a  steady-state 
condition  will  be  obtained.  If  the  stiffness  is  not  adequate  for  the  flow  condition,  the  iterative  process  will 
not  teach  a  steady-state  condition. 
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V.  Discussion  and  Results 


The  results  of  the  research  are  presented  as  solutions  of  potential,  velocity,  coefficient  of  pressure 
and  deformation.  Because  the  solution  of  potential  flow  about  a  cylinder  is  well  known,  the  results  of  the 
BEM  potential  flow  are  compared  to  theory  to  validate  the  potential  flow  code.  Once  the  potential 
solutions  have  been  discussed,  the  results  of  the  elastic  deformations  are  discussed  and  are  compared  to 
results  from  a  FEM  solution.  Finally,  the  iteration  of  the  fluid-structure  problem  is  discussed. 

Effect  of  Element  Size  on  Solution 

The  structure  could  be  modeled  with  8  to  120  elements,  and  as  the  number  of  elements  increased  the 
size  of  each  element  decreased.  As  the  element  size  decreased,  the  geometry  approached  that  of  a  cylinder. 
Therefore,  it  would  be  expected  that  as  the  size  of  elements  decreased  the  BEM  solution  for  the  velocity  and 
pressure  would  approach  theoretical  values  about  a  cylinder. 

The  coefficient  of  pressure  was  the  parameter  chosen  to  discuss  because  it  is  the  end  result  of  the 
potential  and  velocity  calculations.  Figures  Al,  A2  and  A3  illustrate  the  effect  of  element  size  on  the 
solution  of  coefficient  of  pressure  (Cp)  for  constant,  linear  and  parabolic  modeling  of  potential  and 
velocity.  The  solution  is  presented  for  8  elements,  40  elements  and  120  elements  about  the  structure.  The 
angle  around  the  cylinder.  Theta  (8),  is  measured  from  the  front  stagnation  point  in  a  clockwise  direction 
(Figure  14).  The  constant  order  approximation  results  in  step  function  solutions,  the  linear  approximation 
as  a  linear  solution  from  one  end  of  the  element  to  the  other,  and  the  parabolic  approximation  results  in 
parabolic  behavior  of  the  solution.  The  spikes  in  the  parabolic  solution  are  expected  due  to  the 
discontinuity  of  the  geometry  slope  as  discussed  in  Chapter  4  and  in  Reference  8.  Note  that  as  the  element 
size  decreases,  the  amplitude  of  the  spikes  decreases  also,  although  the  ratio  of  element  size  to  spike 
amplitude  increases.  If  the  geometry  was  modeled  with  slope  continuity,  it  would  be  expected  that  the 
parabolic  approximation  solution  spikes  would  smooth. 
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Figure  14.  Definition  of  Angle  Around  Cylinder 


Effect  of  Approximation  Order  on  Solution 

The  approximations  are  compared  against  each  other  for  8, 40  and  120  elements  in  Figures  A3,  A4 
and  A5.  Higher  order  approximations  than  parabolic  would  also  display  singular  behavior  near  the  comers. 

Effect  of  End  Node  Location  on  Solution 

The  effect  of  moving  the  end  nodes  towards  the  element  comer  are  shown  in  Figures  A7  and  A8. 
Parabolic  approximating  order  was  used  for  this  comparison.  The  location  of  the  end  nodes  does  appear  to 
effect  the  solution.  As  shown,  moving  the  nodes  closer  to  the  comers  causes  the  spiking  behavior  to 
amplify.  This  is  due  to  the  ln(l/r)  and  1/r  singular  behavior  of  the  fundamental  solutions  which  cause  the 
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integrals  to  increase  towards  infinity.  Indeed,  when  the  nodes  were  placed  within  .1  5  of  the  comers,  the 
solution  quickly  became  invalid. 


Comparison  of  Potential  Flow  Results  With  Theory 

Linear  approximations  with  120  elements  were  chosen  for  the  comparison  of  the  BEM  results  with 
theory.  Figures  A9  through  A12  show  comparisons  of  the  potential  flow  velocities  and  Cp  calculations. 
For  an  fluid  flow  velocity  of  1,  theoretical  values  of  Cp  are  1  at  the  front  and  back  stagnation  points 
(0  =  0, 180),  and  -3  as  the  flow  acclerates  around  the  top  and  bottom  of  the  cylinder  (0  =  90, 270).  As 
shown,  the  results  of  the  code  model  the  stagnation  pressure  well  with  any  element  size.  The  improvement 
of  the  solution  due  to  decreasing  element  size  is  most  apparent  at  the  top  and  bottom  of  the  cylinder  as  the 
solution  converges  towards  theory.  The  largest  error  of  the  BEM  appears  at  the  top  and  bottom  of  the 
structure  and  is  on  the  order  of  2.4%  of  Cp.  Figure  A13,  which  displays  the  comparison  of  BEM  potential 
with  theory,  shows  that  the  BEM  results  overlay  those  of  theory. 

Effect  of  Number  of  Elements  and  Approximation  Order  on  Solution  Time 

The  elastic  deformations  were  calculated  using  40  flat  elements  and  a  linear  approximation  of  the 
potential  and  deformation.  Although  120  elements  produced  a  slightly  smoother  pressure  distribution  than 
40  elements,  it  also  significantly  increases  the  computer  time  required  for  solution.  This  is  because  a 
(8xm)  x  (8xm)  system  of  equations  is  solved  to  determine  the  deformation.  For  40  elements  with  linear 
modeling  this  produces  a  640  x  640  system  of  equations  which  took  approximately  40  minutes  to  execute 
on  the  AFIT  Interim  Computer  Capability  (ICC)  system.  The  ICC  is  an  Elxsi  6400  superminicomputer 
with  the  Berkeley  4.3  UNIX  operating  system.  A  1920  x  1920  system  of  equations  is  produced  for  120 
elements  with  linear  modeling.  Due  to  this  factor,  it  was  felt  that  the  solution  with  40  elements  was 
adequate  to  show  trends  and  results  of  the  BEM  code.  For  flat  elements  and  linear  approximation  of 
deformation,  a  larger  number  of  elements  would  serve  to  slightly  smooth  the  final  results. 
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Elastic  Deformation  Results 


The  density  of  the  potential  fluid  was  defined  as  that  of  water  at  standard  day  temperature  (1.937 
slugs/ft-3).  The  structure  was  modeled  with  a  Young's  modulus  of  1000,  and  a  Poisson's  ratio  of  0.3. 
These  assignments  were  made  so  that  the  resulting  deformations  would  be  large  and  more  apparent 

The  x  and  y  components  of  deformation  about  the  outer  boundary  are  shown  in  Figures  A 14  and 
A15.  The  figures  represent  the  structure  being  deformed  inward  on  the  front  and  back  sides,  and  pulled 
outward  at  the  top  and  bottom  (Figure  15).  rhis  outward  pulling  is  caused  by  the  accelerated  flow 
producing  negative  values  of  Cp,  and  therefore  the  pressure  at  the  top  and  bottom  is  outward.  The 
deformation  along  the  inner  boundary,  although  not  shown,  was  outward  in  all  directions,  nearly  equally 
opposing  the  outer  boundary  on  the  front  and  back  and  pushing  outward  in  the  same  amount  as  the  outer 
boundary  on  the  top  and  bottom  (Figure  15).  The  inflections  of  the  deformation  curves  are  due  to  the 
element  bulging. 
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Rigid-Bodv  Constraint 

Note  that  although  the  structure  was  not  physically  pinned,  rigid-body  movement  is  constrained  and 
the  structure  does  not  translate  or  rotate.  There  is  zero  deformation  in  the  x  direction  at  the  top  and  bottom 
restricting  side-to-side  movement,  and  zero  defamation  in  the  y  direction  at  the  sides  retricting  up  and  down 
movement.  This  is  due  to  the  method  of  calculating  the  coefficient  c1  with  consideration  of  rigid-body 
movements.  Effectively,  by  calculating  the  diagonal  terms  of  the  H  matrix  from  a  summation  of  the  terms 
on  the  row,  the  solution  is  constrained.  Additional  boundary  conditions  may  be  placed  on  the  structure  if 
desired,  and  the  system  of  equations  treated  accordingly. 

Comparison  With  Finite  Element  Model 

A  finite  element  mesh  was  constructed  with  40  linear  elements.  The  internal  pressure  calculated 
from  the  BEM  potential  flow  code  was  distributed  at  all  the  internal  nodes.  Concentrated  forces  were  placed 
at  the  outer  boundary  nodes  to  represent  the  pressure  distribution  determined  from  the  BEM  solution. 

Because  the  FEM  nodes  were  at  the  comers  of  the  elements,  the  BEM  pressure  was  calculated  at  the 
end  points  using  shape  functions  and  averaged  between  elements  to  determine  the  concentrated  pressure 
values  at  the  comers. 

The  deformations  from  the  FEM  code  are  shown  in  Figures  A16  and  A17.  The  shape  of  the  BEM 
and  FEM  solutions  are  nearly  identical,  and  although  the  magnitudes  are  different,  they  are  consistent  with 
past  experience  of  the  FEM  code  used. 

Iteration  of  Fluid/Structure  Solutions 

A  linking  geometry  routine  was  written  to  apply  the  elastic  deformations  to  the  original  structure 
boundary,  and  then  to  recalculate  the  potential  and  structural  solutions.  Once  the  original  geometry  was 
defined,  no  additional  information  was  required  from  the  user  to  iterate  on  the  solutions.  This  factor  is  the 
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major  advantage  to  a  boundary  element  formulation  over  the  domain  formulations  as  discussed  in  this 
research. 

For  the  iterative  results,  8  elements  were  used  with  linear  approximation  functions.  The  first 
calculated  x  and  y  components  of  deformation  are  plotted  in  Figures  A18  and  A20.  Results  from  the  third 
to  the  9th  iterations  are  shown  in  Figures  A 19  and  A21.  As  shown,  the  magnitudes  of  the  deformations 
decrease  with  each  interation,  and  have  approached  zero  within  12  decimal  places  by  the  9th  iteration 
thereby  indicating  a  steady-state  solution  has  been  obtained. 
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-VI  Conclusions  and  Recommendations 

The  major  objective  of  this  investigation  was  to  determine  the  application  of  the  boundary  element 
technique  to  compute  the  incremental  deformations  of  an  elastic  structure  in  fluid  flow.  The  main 
advantage  to  the  boundary  dement  formulation  over  finite  differences  or  finite  elements  is  that  it  can 
numerically  approximate  both  the  fluid  and  elastic  structure  governing  equations  with  a  common  definition 
of  the  fluid/structure  boundary. 

Results  obtained  during  the  course  of  this  research  indicate  the  the  BEM  does  have  a  favorable 
application  to  the  aeroelasdc  problem  of  fluid  flow  past  an  elastic  body.  While  simplifying  assumptions 
have  been  made  in  the  execution  of  this  thesis,  it  could  be  seen  that  this  method  presents  a  user  friendly 
approach  to  the  problem  definition.  The  fact  that  once  the  original  boundary  is  specified  the  user  is  not 
required  to  continuously  redefine  the  geometry  would  be  a  large  advantage  in  using  organizations  such  as 
the  4950th  Test  Wing.  Also,  as  implemented  in  this  study,  the  user  is  able  to  use  several  different 
approximating  orders  rather  than  be  limited  to  one.  This  would  allow  better  modeling  of  the  solution  as  a 
function  of  the  problem  to  be  solved. 

As  shown  with  the  results  from  the  program,  the  BEM  can  provide  accurate  results.  The  iterative 
capability  of  the  code  was  also  verified.  This  capability  should  prove  to  be  an  improvement  over  current 
computational  techniques  of  the  calculation  of  airloads  and  deformations  in  an  aeroelastic  system.  While 
the  process  did  produce  steady-state  solutions,  the  stability  of  the  BEM  iterative  formulation  should  be 
more  thoroughly  researched.  Stability  criteria,  such  as  that  used  in  finite  difference  and  finite  element 
codes,  may  become  applicable. 

The  largest  factor  towards  accurate  implementation  of  the  formulation  is  the  approximation  of  the 
boundary  geometry.  While  presenting  a  simple  method  of  implementation,  discontinuous  geometry  slope 
produces  many  undesirable  results  and  the  use  of  continuous  higher  order  geometry  elements  is  highly 
recommended  for  future  researchers. 

Several  additional  features  could  be  added  to  increase  the  versitility  of  the  code.  For  example; 
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three-dimensional  problems,  adding  wakes  to  properly  solve  for  lifting  body  problems,  and  using 
non-stead;  problem  implementation.  The  applicability  of  the  FEM  to  aeroelastic  problems  would  be 
increased  by  each  of  these  additions. 
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Appendix  A:  Investigation  Results 
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Figure  Al.  Effect  of  Element  Size  on  Solution,  Constant  Modeling  of  Cp 
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Figure  A3.  Effect  of  Element  Size  on  Solution,  Parabolic  Modeling  of  Cp 
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Figure  A4.  Effect  of  Function  Approximation  on  Solution,  8  Elements 
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Figure  A5.  Effect  of  Function  Approximation  on  Solution,  40  Elements 
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Figure  A6.  Effect  of  Function  Approximation  on  Solution,  120  Elements 


Figure  A7.  Effect  of  End  Node  Location  on  Solution,  40  Elements,  Parabolic  Approximation  of  Cp 
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Figure  A8.  Effect 
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Figure  A9.  Comparison  of  X  Component  of  Velocity  With  Theory,  120  Elements,  Linear  Approximation 


Comparison  of  Y  Component  of  Velocity  With  Theory,  120  Elements,  Linear  Approximation 


Figure  A12.  Comparison  of  Coefficient  of  Pressure  With  Theory,  120  Elements,  Linear  Approximation 
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Figure  A13.  Comparison  of  Potential  With  Theory,  120  Element,  Linear  Approximation 
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Figure  A 14.  X  Component  of  Deformation,  40  Elements,  Linear  Approximation 
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Figure  A 16.  X  Component  of  Deformation  from  FEM,  40  Elements,  Linear  Approximation 
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Figure  A 17.  Y  Component  of  Deformation  From  FEM,  40  Elements,  Linear  Approximation 


Figure  A18.  X  Component  of  Deformation  After  1  Pass,  8  Elements,  Linear  Approximation 
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Figure  A 19.  Iterative  X  Component  of  Deformation,  3-9  Iterations,  8  Elements,  Linear  Approximation 
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Figure  A20.  Y  Component  of  Deformation  After  1  Pass,  8  Elements,  Linear  Approximation 


Appendix  B:  Program  Routines 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

C  FLUID  -  STRUCTURE  INTERACTION  USING  BOUNDARY  ELEMENT  TECHNIQUE 
C 

C  THESIS  PROJECT 
C  NORMA  F  TAYLOR 
C  GAE-88D 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 

c 

C  MAIN  PROGRAM 
C 

C  IF LAG  -  1  PERFORM  FLUID  CALCULATIONS  ONLY 
C  -  OTHER  ITERATE  BETWEEN  FLUID  AND  STRUCTURE 

C  PFLAG  -  1  PRINT  OUT  POTENTIAL  INFORMATION 

C  -  2  PRINT  OUT  STRUCTURAL  INFORMATION 

C  -  3  PRINT  OUT  BOTH 

C 

implicit  real*8  (a-h,o-t) 

COMMON /FLAGS/IFLAG, PFLAG, ITER 

COMMON /GEOMET/NNO DES , NELEMS, ORDER, EPS , NORM I , NORMJ , JAC, NODEC , 
•NODPIN 1 , NODPIN  2 ,  NODE ( 5 ) , ELEM (360,4) 

COMMON/POTEN/VINF , ALPHA , VI , VJ , PINF , RHOINF , GAMMA , MINF2 , NODEP 
COMMON/FLOW/H (240 , 240 ) ,G ( 240, 240 ) ,Q ( 240 ) , GHAT ( 240  ) 
COMMON/STRUCT/GE,XNU,HS( 1920, 1920 ) ,GS( 1920, 1440 ) ,GSHAT( 1920  ) 
COMMON /POSITION /XI ( 1080 ) , YI ( 1080 ) 

COMMON /VELOCIT/UVEH  240  ) ,  WEL{240  )  ,XFORC ( 480  ) , YFORC  ( 480  ) , 

•XOLD ( 480 ) , YOLD ( 480 ) 

REAL*8  NODE, GAMMA, MINF 2, NORM I, NORMJ, JAC 
INTEGER  ORDER, PFLAG 
C 
C 
C 

R£AD(5, 1 )  IFLAG, PFLAG 

1  FORMAT (211) 

READ( 5, 2 )  ITERATE 

2  FORMAT (12) 

C 

IF( IFLAG. EQ. 1 )  THEN 
C 

C  PERFORM  FLUID  CALCULATIONS  ONLY 
C 

C  READ  IN  GEOMETRY  AND  FLUID  INFORMATION 
C  FORM  H ( PHI )-GHAT  SYSTEM 
C  SOLVE  FOR  PHI 

CALCULATE  PRESSURES  (XFORC, YFORC  ) 

CALL  GEOMIN 
CALL  GOUTPUT ( ITER ) 

CALL  FLU I DIN 
CALL  FFORM 
CALL  QFLOW 
DO  5  I" 1, NNODES 
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no  onnnnnooonooo 


MAIN  2/3 


GHAT(I)-0.0 
DO  5  J- 1 , NNODES 

GHAT(I )K3HAT(I )+G(I,J)*Q(J) 

5  CONTINUE 

CALL  HDIAG 

IF(PFLAG.EQ.  1.0R.PFLAG.EQ.3)  CALL  HGQOOT 
CALL  SOLVE <G,GHAT,D, NNODES, 240  ) 

C 

C  NOTE  THAT  AFTER  GOING  THROUGH  SOLVE  GHAT  IS  THE  SOLUTION  PHI 
C 

CALL  VELOC 
CALL  PRESSURE 
CALL  FLOWOOT 
CALL  DRAW (ITERATE) 

ELSE 

CALCULATE  BOTH  FLUID  AND  STRUCTURAL  INFORMATION 
AND  ITERATE  UNTIL  STEADY  STATE  SOLUTION  IS  REACHED 

READ  IN  GEOMETRY,  FLUID  AND  STRUCTURAL  INFORMATION 
CALCULATE  H( PHI) -GHAT  SYSTEM 
SOLVE  FOR  PHI 

CALCULATE  SURFACE  PRESSURES 
CALCULATE  HS(U)-GSHAT  SYSTEM 
SOLVE  FOR  DEFORMATIONS  U 
CALCULATE  NEW  GEOMETRY 
ITERATE 

ITER-0 

CALL  GECMIN 

CALL  GOUTPUT ( ITER ) 

CALL  FLUIDIN 
CALL  SINFO 
100  CALL  FFORM 
CALL  QFLOW 
DO  105  1-1, NNODES 
GHAT (I )-0.0 
DO  105  J-1, NNODES 

GHAT (I )-GHAT(I )+G  ( I ,  J )  *Q  ( J  ) 

105  CONTINUE 

CALL  HDIAG 

IF(PFLAG.EQ. 1 . OR. PFLAG. EQ. 3 )  CALL  HGQOOT 
CALL  SOLVE(H, GHAT, D, NNODES, 240) 

NOTE  THAT  AFTER  GOING  THROUGH  SOLVE  GHAT  IS  REALLY  THE  SOLUTION  PHI 

CALL  VELOC 
CALL  PRESSURE 
CALL  FLOWOUT 
CALL  SFORM 
CALL  HSDIAG 
CALL  GSKATT 

IF ( PFLAG, EQ, 2 . OR. PFLAG . EQ. 3 )  CALL  SOOT 
CALL  SOLVE ( HS, GSHAT, D, NNODES *8, 1920 ) 
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MAIN  3/3 
C 

C  NOTE  THAT  AFTER  SOLVES  GSHAT  IS  REALLY  THE  DEFORMATION  U 
C 

CALL  SOUT PUT ( ITER ) 

ITER-ITER+1 
CALL  DRAW (ITER) 

CALL  DEFENDS 
CALL  GOUTPUT(ITER) 

IF( ITER. LT. ITERATE)  GO  TO  100 
ENDIF 

150  CONTINUE 
END 
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ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

C  DEFENDS 

c 

C  THIS  ROUTINE  CALCULATES  THE  DEFORMATIONS  AT  THE  ENDS  OF  THE  ELEMENTS 
C  AND  DETERMINES  THE  NEW  ELEMENT  COORDINATES 
C 

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC 
SUBROUTINE  DEFENDS 
IMPLICIT  REAL*8  (A-H.O-Z) 

COMMON /GEOMET/NNODES , NELEMS , ORDER , EPS , NORM I , NORMJ , JAC , NODEC , 

'NODPIN 1 , NOD PIN 2 , NODE ( 5 > , ELEM ( 36 0 , 4 ) 

COMMON/STRUCT/GE , XNU, HS (1920,  1920  )  ,GS ( 1920, 1440),GSHAT( 1920) 

COMMON /POSIT  ION /XI ( 1080 ) ,YI ( 1080 ) 

REALMS  NODE,NORMI,NORMJ,JAC,SHAP(5),SHP(5) 

INTEGER  ORDER 
C 

C  CALCULATE  THE  X  AND  Y  DEFORMATION  OF  THE  LEFT  END  OF  ELEMENT  J 
C  BY  AVERAGING  THE  DEFORMATION  AT  ETA— 1  E’ EMENT  J  AND 
C  ETA-1  ELEMENT  J-1 
C 

DO  20  J— 1 , 2*NELEMS 
UXL-0.0 
UXR-0 . 0 
UYL-0.0 
UYR-0.0 

CALL  SHAPE (ORDER, -1. DO, SHAP) 

CALL  SHAPE(ORDER, 1.D0,SHP) 

DO  10  L-1, ORDER 

JL-J 'ORDER- (ORDER-L ) 

JR- (J-1 ) 'ORDER- (ORDER-L) 

IF(J.SQ.I)  JR-NELEMS 'ORDER- (ORDER-L ) 

IF( J.EQ.NELEMS+1 )  JR— 2*NELEMS*ORDER-(ORDER-L ) 

UXL-UXL+SHAP  ( L )  'GSHAT  ( JL ) 

UXR-UXR+S  HP ( L ) 'GSHAT (JR) 

UYL-UYL+SHAP (L )  'GSHAT ( 4 'NNODES+ JL ) 

UYR— UYR+SHP  ( L )  *GS  HAT  ( 4  *NNODES  +JR ) 

10  CONTINUE 

ELEM( J , 1 )-ELEM( J , 1 )  +  (UXL-KJXR  )/2 . 

ELEM( J, 2 )-ELEM( J, 2 )+ ( UYL+UYR ) /2 . 

20  CONTINUE 
C 

C  CALCULATE  THE  RIGHT  ENDS  OF  THE  ELEMENTS  BASED  ON  THEM 
C  EQUALING  THE  LEFT  ENDS  OF  THE  J+1  ELEMENT 
C 

DO  30  J-1, 2 'NELEMS 
IF(J.BQ. NELEMS)  THEN 
ELEM ( J, 3) -ELEM ( 1,1) 

ELEM (J, 4 )-ELEM( 1,2) 

ELSE 

IF(J.EQ.2*NELEMS  )  THEN 

ELEM( J, 3 ) -ELEM (NELEMS -M, 1  ) 

ELEM< J, 4 )-ELEM( NELEMS +  1 , 2  ) 
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DEFENDS  2/2 


ELSE 

ELEM{ J , 3 )-ELEM(J+1 , 1 ) 

ELEM(  J,4  )-ELEM(J-H,2  ) 

END  IF 
END  IF 
CONTINUE 

CALCULATE  THE  CONNECTING  ELEMENT  END  POINTS 

DO  40  I“1 , NELEMS 

ELEM(2*NELEMS+I , 1 )-ELEM< I ,3) 

ELEM(2 ‘NELEMS+I , 2 ) -ELEM (1,4) 

ELEM( 2*NELEMS+I , 3 )«ELEM(NELEMS+I , 3  ) 

ELEM(2*N£LaMS+I,4 )«ELEM( NELEMS+I , 4 ) 

CONTINUE 

CALCULATE  THE  NEW  COORDINATES  OF  THE  NODAL  POINTS 

DO  50  I*1,NELEMS*3 
DO  50  L«1, ORDER 
ETAA— 1  . 

ETAC-1. 

ETAQ-NODE(L) 

SHAP 1 - ( ETAQ-ETAC ) / ( ETAA-ETAC ) 

SHAP2- ( ETAQ-ETAA ) / ( ETAC-ETAA ) 

XI ( I ‘ORDER- (ORDER-L ) ) *SHAP 1 *ELEM ( I , 1 )+SHAP2  *ELEM (1,3) 
YI ( I ‘ORDER- ( ORDER-L ) ) *SHAP 1 ‘ELEM (1,2) +SHAP2  *ELEM (1,4) 
CONTINUE 
RETURN 
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.............. ••••••••• 

C 

C  FFORM 
C 

C  THIS  SUBROUTINE  COMPUTES  THE  H  AND  G  MATRICES  FOR  THE 
C  POTENTIAL  FLOW 
C 

C  H  MATRIX  -  INTEGRAL(SHAPE*QSTAR*JAC)DETA 
C 

C  G  MATRIX  -  INTEGRAL (SHAPE *PH ISTAR*JAC )DETA 
C 

c.. 

SUBROUTINE  FFORM 
implicit  real'8  (a-h,o-z) 

COMMON/POSITION/XI { 1 080 ) , YI ( 1080) 

COMMON /GEOMET/NNODES , NELEMS , ORDER, EPS , NORM I , NORMJ , JAC , NODEC , 
'NODPIN 1 , NODPIN2, NODE ( 5 ) , ELEM( 360, 4 ) 

COMMON /FLOW/H< 240, 240 ), G< 240, 240 ),Q( 240), GHAT (240) 
real *8  ETA ( 12 ) , WT ( 12), ETALN  < 10 ) , WTLN ( 10 ) , ETA2 ( 10 ) ,WT2 ( 10 ) 
real*8  SHAP( 5 ) ,SHPLN2( 5 ) , SHPLN(5 ) ,SHAP3 ( 5 ) 
real'8  NODE, JAC, JACLN, JACLN2, JAC3, JAC4,NORMI , NORMJ 
INTEGER  ORDER 
C 

C  ETA  AND  WT  ARE  NORMAL  12TH  ORDER  GAUSSIAN  QUADRATURE 

C  VALUES  AND  WEIGHTS 

C 

DATA  ETA( 1 ) , ETA < 2 ) ,ETA( 3 ) ,ETA( 4 ) , ETA ( 5 ) , ETA ( 6 ) , ETA ( 7  )  , 

'ETA ( 8 ) , ETA ( 9 ) , ETA ( 10 ) , ETA< 1 1 ) , ETA ( 12 ) /-. 981 56063d0 , -. 9041 1 726d0 , 
76990267d0 , 5873 1 795d0 ,-. 36783 149d0 , -.  1 252  334 1 dO , 

*. 12523341 dO, . 36783 149d0 ,. 5873 1795d0 , . 76990267d0, .9041 1726d0, 

* .98 1 56063d0/ 

DATA  WT ( 1 ) , WT{2 ) , WT( 3 ) , WT( 4 ) , WT( 5 ) , WT(6 ) ,WT( 7 ) , WT( 8 ) , 

*WT(9),WT( 10),WT( 11,,WT( 1 2 )/. 0471 7534d0, . 10693933d0 , . 16007833d0 , 
*.20316743d0, .23349254d0, .24914705d0, .24914705d0, .23349254d0, 

*. 20316743 d0,.1 6007 83 3 dO, . 10693933d0, . 04717534d0/ 

C 

C  ETALN  AND  WTLN  ARE  10TH  ORDER  LOGARITHMIC  GAUSSIAN  QUADRATURE 

C  VALUES  AND  WEIGHTS 

C 

DATA  ETALN ( 1 ) , ETALN ( 2 ) , ETALN ( 3 ) , ETALN  <  4 ) , ETALN ( 5 ) , ETALN ( 6  ) , 
'BTAIN  ( 7 ) , ETALN ( 8 ) , ETALN ( 9 ) , ETALN (10)/. 009042 59 dO, .0  5397105d0, 

*. 135311 34d0, . 2470 5 169d0, .38021 171 dO, . 523791 59d0, .66577472d0, 

' . 794 1 90 19d0, .89816102d0, .96884798d0/ 

DATA  WTLN ( 1 ) , WTLN ( 2 ) , WTLN ( 3 ) , WTLN (4 ) , WTLN ( 5 ) , WTLN (6 ) , WTLN ( 7 ) , 
•WTLN (8 ) ,WTLN ( 9 ) , WTLN  (  10  )/.  120  95474d0, .  1 863631 OdO, . 19566066d0, 

'. 17357723d0, . 13569597d0, . 09364 7C84dO , . 055 787938d0 , . 027 1 59893d0 , 
*. 00951 51 992 dO, .00 16381 586d0/ 

C 

C  ETA 2  AND  WT2  ARE  10TH  ORDER  GAUSSIAN  QUADRATURE 
C  VALUES  AND  WEIGHTS 
C 

DATA  ETA2 ( 1 ) ,ETA2 ( 2 ) ,ETA2 ( 3  )  , ETA 2 ( 4 ) , ETA 2(5), ETA 2(6), ETA 2(7), 
*KTA2 ( 8 ) , ETA2 ( 9 ) , ETA 2 (10)/-.97390653d0,-.86506337d0,-.67940957d0, 
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o  n 


FFORM  2/4 

*-.43339539d0, 14887434d0, .  14887434d0, .43339539d0, .67940957d0, 
*.86506337d0,.97390653d0/ 

DATA  WT2(1 ).WT2(2),WT2<3},WT2(4),WT2(5),WT2(6),WT2(7),WT2{8), 
*WT2(9),WT2( 10 )/. 06667 134d0, . 149451 35<J0 , . 21 908636d0 , . 26926672d0 , 
*.29552422d0, .29552422d0, .26926672d0, .21908636d0, . 14945135d0, 
*.06667134d0/ 

DATA  PI/3. 14159265d0/ 

C 

C  INITIALIZE  MATRICES 
C 

DO  10  1*1 , NNODES 
DO  10  J-1.NN0D8S 
G(I,J)-0.0 
H(I,J)-0.0 
10  CONTINUE 
C 

C  BEGIN  LOOPS  ON  NODE  I  AND  ELEMENT  J 

C  FOR  EACH  I  ROW,  PERFORM  INTEGRATIONS  OVER  EACH  ELEMENT  J 
C  AROUND  THE  STRUCTURE.  FILL  IN  ONE  ROW  AT  A  TIME 
C 

DO  1000  1-1, NNODES 
DO  1000  J- 1 , NELEMS 
C 

C  DETERMINE  IF  THE  INTEGRATION  SHOULD  BE  SINGULAR 
C  I.E.  IS  THE  NODE  I  SOMEWHERE  ON  THE  ELEMENT? 

C  ISING  -  0  MEANS  NODE  I  IS  NOT  ON  ELEMENT  J 

C  ISING  -  1  MEANS  NODE  I  IS  FIRST  NODE  ON  ELEMENT  J 

C  ISING  -  2  MEANS  NODE  I  IS  SECOND  NODE  ON  EliMENT  J 

C  ISING  -  3  ETC 

C 

ISING— 0 

DO  20  L-1, ORDER 

NODES-J 'ORDER- (ORDER-L ) 

IF  ( I  ■> EQ. NODES )  ISING-L 
20  CONTINUE 

IF( ISING. EQ. 0 )  THEN 

C  THE  INTEGRATION  IS  NOT  SINGULAR 
C  PERFORM  12TH  ORDER  GAUSSIAN  QUADRATURE 

•»*.***••  •**•••***»**»•»* 

DO  50  K-1,  12 

CALCULATE  THE  »  HD  Y  COORDINATES  OF  THE  QUADRATURE  POINT 
C  AND  CALCULATE  DISTANCE  FROM  NODE  I  TO  QUAD  POINT 
C 

C  CALCULATE  SHAPE  FUNCTIONS  AT  QUAD  POI..T  ( ETA  ( K )  )  BASED  ON  GEOMETRY 
C  NODES  AT  -1,0,1 
C 

ETA  A—  1 . 

ETAC-1. 

ETAQ-ETA ( K ) 

SHAP( 1 )- ( ETAQ-ETAC ) /( ETAA-ETAC ) 

S HAP  (  2  ) - ( ETAQ -ETAA  )  /( ETAC -ETAA  ) 

X-SHAP( 1 )*ELEM(J, 1 )+SHAP( 2 ) *ELEM( J , 3 ) 
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Y -SHAP (1 ) *ELEM ( J , 2 ) +SHAP ( 2 ) *ELEM ( J , 4 ) 

RI-XI (I )-X 
RJ-YKD-Y 
R«SQRT(RI*RI+RJ*RJ) 

C 

C  CALCULATE  PHI*  AND  Q* 

C 

CALL  JACOB IAN (J,ETA(K) ) 

PHISTAR-dlogO/R) 

D-RI * ( NORM I ) +RJ  * ( NORM J ) 

QSTAR-(D/(R*R) ) 

C 

C  CALCULATE  TERMS  OF  H  AND  G  MATRICES 
C  ROW  CORRESPONDS  TO  NODE  I 

C  COLUMN  CORRESPONDS  TO  GLOBAL  NUMBER  OF  ELEMENT  NODES 
C 

CALL  SHAPE ( ORDER, ETA(K) , SHAP) 

DO  35  L-1, ORDER 
C 

C  JJ  IS  GLOBAL  NUMBER  OF  ELEMENT  NODES 
C  FOR  EXAMPLE;  NODES  1,2,3  ON  ELEMENT  1  ARE  1,2,3 

C  NODES  1,2,3  ON  ELEMENT  2  ARE  4,5,6 

C  NODES  1,2,3  ON  ELEMENT  3  ARE  7,8,9 

C 

J J“J *0 RDER “( ORDER -L ) 

G(I,JJ)-G(I,JJ)+SHAP(L)*PHISTAR*JAC*WT(1C)/<2.*PI) 
H { I , JJ ) «H ( I , J J ) +SHAP < L ) *QSTAR* J AC*WT ( K ) / ( 2 . *PI ) 

35  CONTINUE 

50  CONTINUE 

ELSE 

c 

C  THE  INTEGRAL  IS  SINGULAR  IN  THE  INSIDES  OF  ELEMENT  J 
C 

C  DETERMINE  THE  ETA  BAR  LOCATION  OF  SINGULARITY 
C 

8 INGUL-NODE ( I S ING ) 

C 

C  CALCULATE  TERMS  OF  G  MATRIX 

C 

C 

X2X 1 “ELEM( J , 3 )-ELEM( J , 1 ) 

Y2Y 1 -ELEM( J, 4 ) -ELEM( J, 2 ) 

PLENG«SQRT(X2X1*X2X1+Y2Y 1*Y2Y 1 ) 

DO  350  X-1,10 

ET-S INGUL+ ( 1 . -S INGUL ) *ETALN ( K ) 

CALL  SHAPE (ORDER, ET, SHPLN2 ) 

CALL  JACOBIAN (J,ET) 

JACLN2-JAC 

ET-SINGUL+.5*( 1 . -S INGUL )•( ETA2 (K ) -M  . ) 

CALL  SHAPE ( ORDER, ET, SHAP ) 

CALL  JACOBIAN) J, ET  ) 

JAC4-JAC 

BT-SINGUL-(S INGUL*  1.  )*ETALN(JC) 
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CALL  SHAPE ( ORDER, ET, SHPLN ) 

CALL  JACOB IAN (J,ET) 

JACLN-JAC 

ET-SINGUL-.5*( 1.+SINGUL)*(ETA2(K)+1. ) 

CALL  SHAPE (ORDER, BT, SHAP3 ) 

CALL  JACOBIAN <J,ET) 

J  AC  3  “J  AC 

DO  385  L“1, ORDER 

JJ-J ‘ORDER- ( ORDER-L ) 

TERM  1 “SHPLN (L ) *JACLN*( SINGUL+1 . ) *WTLN (K ) /( 2 . *PI ) 
TERM3“SHAP3 (L) *JAC3* 

dlog(PLENG*ABS(1.+SINGUL)/2. )*(SINGUL+1. ) 
*.5*WT2(K)/(2.*PI) 

TERM2«SHPLN2(L)*JACLN2*( 1 . -SINGUL ) *WTLN (K )/( 2 . *PI ) 
TERM4-SHAP(L)*JAC4* 

dlog( PLENG*ABS ( SINGUL- 1 . )/2. )*( 1 .-SINGUL) 

* . 5*WT2  <K )/ ( 2 . *PI ) 

G(I,JJ)<(I,  JJ  )+TERM 1 +TERM2-TERM3-TERM4 
CONTINUE 
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350  CONTINUE 

END  IF 

1000  CONTINUE 
RETURN 
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C 

C  FLOWOUT 
C 

C  THIS  SUBROUTINE  OUTPUTS  THE  POTENTIAL  FLUID  INFORMATOINS 
C 

C***************************************************************** 

SUBROUTINE  FLOWOUT 
implicit  re«l*8  (a-h,o-z) 

COMMON /GEOMET/NNODES , NELEMS , ORDER , EPS , NORMI , NORMJ , JAC , NODEC , 
•NODPIN 1 , NODPIN2, NODE ( 5 ) , ELEM (360,4) 

COMMON/FLOW/H(240,240) ,G(240,240) ,Q(240 ) ,GHAT(240) 

COMMON /POTEN /V IN  F , ALPHA , VI, V J , P I NF , RHO INF , GAMMA , M INF  2 , NODEP 
COMMON/VELOCIT/UVEH240)  ,WEL(240)  ,XFORC(480  )  ,YFORC(480  ), 
•XOLD ( 480 ) , YOLD ( 480 ) 

REAL*8  NORMI , NORMJ , JAC , NODE , GAMMA, M INF 2 
INTEGER  ORDER 


WRITE ( 6 , 1 )  VINF,ALPHA,VI,VJ,MINF2**.5 

1  FORMAT (//, 2 2X, '**•**  FLOW  INFORMATION  *****',/, 24X, 

*  **•  ON  OUTER  BOUNDARY  **',/,  1 0X , 'Vinfinity  *  *,F10.4, 

*  SX, 'ALPHA  -  ' ,F10.4,/, 5X, *VI  -  ' ,F 1 0 . 4 , 1 0X, 'VJ  -  ',F10.4, 

*  10X, 'MINF  -  ' ,  F  1 0 . 4 , // ) 

WRITE (6,2) 

2  FORMAT ( IX, 'POTENTIAL' , 3X, *U  VEL',5X,'V  VEL ',  4X , 'VELOCITY ' , 
•5X, 'CP' ,6X, 'X  FORCE', 3X,'Y  FORCE',//) 

DO  5  I  — 1 , NNODES 

VEL-SQRT(UVEL(I  )**2+WEL(I  )**2  ) 

CP- 1 . - ( VEL/V INF ) * *2 

WRITE (6, 3)  GHAT ( I  )  ,UVEL(I  )  ,WEL(I  )  ,VEL,CP,XFORC(  I  )  ,YFORC(I  ) 

3  FORMAT ( 7F 1 0 . 4 ) 

5  CONTINUE 
C 

C  WRITE  FORCE  ON  INSIDES 
C 

VfiLVENT-SQRT(UVEL(NODEP)**2+WEL(NODEP)**2  ) 

CPVENT- 1 . - ( VELVENT/VINF ) * *2 
PVENT-. 5*RHOINF*VINF*VINF*CPVENT 
WRITE ( 6,f  ;  PVENT, NODEP 

6  FORMAT ( //,22X, ******  FORCES  ON  INNER  SURFACE 

*/, 10X, 'INNER  PRESSURE  IS  -  *,F10.4,/, 

•10X, 'VENTED  NODE  IS  -  ',13,//) 

WRITE (6,7) 

7  FORMAT ( 2X, • X  FORCE', 3X,'Y  FORCE',//) 

WRITE (6, 8)  ( XFORC (I),YFORC(I),I -NNODES +  1 , 2* NNODES ) 

8  FORMAT ( 2  F 1 0 , 4 ) 

RETURN 

END 
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FLU I DIN 

THIS  SUBROUTINE  READS  IN  THE  INFORMATION  REQUIRED  TO 
CALCULATE  THE  POTENTIAL  FLUID  DATA 

INPUT: 

VINF  -  FLOW  VELOCITY  AT  INFINITY 
ALPHA  -  FLOW  ANGLE  OF  ATTACK 
PINF  -  PRESSURE  AT  INFINITY 
RHOINF  -  DENSITY  AT  INFINITY 
GAMMA  -  GAMMA 

NODEP  «  NODE  AT  WHICH  TO  VENT  PRESSURE  TO  INSIDE 
OUTPUT: 

VI  -  X-COMPONENT  OF  FLOW  VELOCITY 
VJ  »  Y-COMPONENT  OF  FLOW  VELOCITY 
MINF2  -  MACH  AT  INFINITY  SQUARED 


SUBROUTINE  FLUIDIN 
implicit  real's  (a-h,o-z) 

COMMON/POTEN/VINF , ALPHA , VI , VJ , PINF , RHOINF , GAMMA, MINF2 , NODEP 
real's  GAMMA, MINF2 

READ  IN  Vinfinity  AND  ANGLE  OF  ATTACK 
CONVERT  Vinf  TO  X  AND  Y  (I  AND  J)  COMPONENTS 

READ( 5 , * )  VINF, ALPHA 
FORMAT (2F1 0.4) 

ALF-ALPHA'O . 0 1 7453293 
VI-VINF'COS(ALF) 

VJ-VINF'SIN(ALF) 

READ  IN  Pinfinity,  RHOinfinity  and  GAMMA 

R£AD( 5, * )  PINF, RHOINF, GAMMA 
FORMAT ( 3F 1 0 . 4 ) 

MINF2«(VINF*VINF*RHOINF)/(GAMMA*PINF) 

READ  IN  NODE  AT  WHICH  TO  VENT  PRESSURE  TO  INSIDE 

READ( 5 , * )  NODEP 
FORMATU2  ) 

RETURN 

END 
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C  GEOMIN 
C 

C  THIS  SUBROUTINE  READS  IN  THE  GEOMETRY  INFORMATION 
C  OF  THE  STRUCTURE 
C 

C  INPUT: 

C  ORDER  -  ORDER  OF  POTENTIAL  AND  Q  MODELING 

C  NELEMS  -  NUMBER  OF  ELEMENTS  ALONG  OUTSIDE  OF  STRUCTURE 

C  TOTAL  NUMBER  OF  ELEMENTS  IS  2*NELEM  BECAUSE  SAME  NUMBER 

C  ALONG  INSIDE  OF  STRUCTURE 

C  EPS  -  DISTANCE  FROM  END  OF  ELEMENT  TO  THE  END  NODES 

C  NODEC  -  NODE  TO  BE  REPLACED  BY  CONTINUITY 

C  NODPIN 1 , NODPIN 2  -  NODES  TO  BE  PINNED 

C 

C  OUTPUT : 

C  ORDER  -  NUMBER  OF  NODES  ON  EACH  ELEMENT 

C  NNODES  -  NUMBER  OF  NODES 

C  NODE ( I )  -  ETA  POSITION  OF  NODE  I  WITHIN  EACH  ELEMENT 

C  XI, YI  -  X,Y  LOCATION  OF  NODE  I 

C 

SUBROUTINE  GEOMIN 
implicit  real*8  (a-h,o-z) 

COMMON/GEOMET /NNODES , NELEMS , ORDER , EPS , NORM I , NORM J , JAC , NODEC , 
•NODPIN 1,NODPIN2,NODE( 5) ,ELEM( 360,4) 

COMMON /POSITION/XI ( 1080 ) ,YI ( 1080  ) 
real *8  NORMI,NORMJ, JAC, NODE 
INTEGER  ORDER 
DATA  PI/3. 14159265d0/ 

C 

READ  IN  ORDER  OF  FIT  AND  NUMBER  OF  ELEMENTS  ALONG  OUTSIDE  SURFACE 

READ( 5 , * )  IORDER, NELEMS 
ORDER-IORDER+1 
NNODES -ORDER*NELEMS 
C 

C  CALL  SURFACE  GEOMETRY  GENERATION  TO  GET  LEFT  AND  RIGHT 
C  X, Y  LOCATIONS  OF  OUTSIDE  AND  INSIDE  PANELS 
C 

CALL  SURF 
C 

C  DETERMINE  LEFT  AND  RIGHT  X,Y  OF  LINKING  PANELS 
C  LINK  OUTSIDE  BOUNDARY  WITH  INSIDE 
C 

DO  10  1-1, NELEMS 

ELEM ( 2  *NELEMS  +  1 , 1 )-ELEM(I,3) 

ELEM ( 2*NELEMS+I ,2 )-ELEM(I,4) 

ELEM( 2*NELEMS+I , 3 )-ELEM( NELEMS+1 , 3  ) 

ELEM ( 2*NELEMS+I , 4 ) =E LEM ( NELEMS + I , 4 ) 

10  CONTINUE 
C 
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READ  IN  DISTANCE  OF  END  NODES  TO  ELEMENT  CORNERS 

READ<5,*)  EPS 
IF(ORDER.EQ. 1 )  EPS-1. 

C  CALCULATE  THE  ETA  POSITION  OF  THE  NODES  WITHIN  THE  ELEMENT 
C  NODE ( 1 )  AND  NODE ( ORDER )  ARE  EPS  FROM  ♦/-  1.  OTHER  NODES 
C  ARE  SPACED  EQUALLY  WITHIN  THE  REMAINING  SPACE 
C 

IF(ORDER. EQ. 1 )  THEN 
NODE ( 1 )-0. 

ELSE 

SPACE— ( 2.-2. *EPS )/( ORDER-1 . ) 

NODE  ( 1  )  — (  1 .  -EPS ) 

NODE  (ORDER)— NODE  (  1  ) 

IF(ORDER.GT.2 )  THEN 
DO  20  1-2, ORDER-1 

NODE ( I ) -NODE ( I - 1 ) +S PACE 
20  CONTINUE 

END  IF 
END  IF 
C 

C  CALCULATE  THE  X.Y  POSITION  OF  EACH  NODE 
C 

DO  40  1—1 ,NELEMS*3 
DO  40  L-1, ORDER 
ETA  A—  1 . 

ETAC-1. 

ETAQ— NODE ( L ) 

SHAP 1 -( ETAQ-ETAC ) / ( ETAA-ETAC ) 

SHAP2- ( ETAQ-ETAA ) / ( ETAC-ETAA ) 

XI ( I *ORDER- (ORDER-L ) J-SHAP1 *ELEM( I , 1 )+SHAP2*ELEM( 1,3) 

Y I ( I ‘ORDER- ( ORDER-L ) ) -S  HAP  1 *ELEM ( 1 , 2 ) +SHAP  2  *ELEM ( 1 , 4 ) 

40  CONTINUE 

C 

C  READ  IN  THE  NODE  OF  EACH  ELEMENT  TO  BE  REPLACED  BY  THE 
C  CONTINUITY  EQUATION 
C 

READ (5,*)  NODEC 
C 

C  READ  IN  NODES  TO  BE  PINNED 
C 

READ (5,*)  NODPIN  1 , NODPIN2 

RETURN 

END 
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C  GOUT PUT 
C 

C  THIS  SUBROUTINE  OUTPUTS  THE  ELEMENT  GEOMETRY  COORDINATES 
C 

c. . . 

SUBROUTINE  GOUT  PITT  (  ITER) 
implicit  real*8  (a-h,o-z) 

COMMON/GEOMET/NNODES , NELEMS , ORDER,  EPS , NORMI , NORMJ , JAC , NODEC , 
•NODPIN 1, NODPIN 2, NODE ( 5 ) , ELEM( 360,4 ) 

COMMON/POSITION /XI ( 1 080 ) , YI ( 1 080 ) 
real*8  NORMI, NORMJ, JAC, NODE 
INTEGER  ORDER 


WRITE (6,10)  NELEMS , 3*NELEMS , ORDER, NNODES , 3 *NNODES , EPS , ITER- 1 
10  FORMAT(//,22X, ***** ‘GEOMETRY  INFORMATION*****1,//, 

*1  OX, ’NUMBER  OF  ELEMENTS  ON  OUTSIDE  -  ’,13, 

*  1 0X, ‘TOTAL  NUMBER  OF  ELEMENTS  -  ',13, 

*/, 10X, ‘ORDER  OF  MODELING  -  *,I3, 

*/, 1 0X, ’NUMBER  OF  NODES  ALONG  OUTSIDE  -  ’,13, 

*10X, ’TOTAL  NUMBER  OF  NODES  -  ’,13, 

*/,1 OX, ’DISTANCE  FROM  END  NODES’, 

*’  TO  CORNER  -  ’ ,F10.4,/, 10X, ’  ITERATION  NUMBER  «  *,I2,/) 
WRITE(6,20) 

20  FORMAT (//,5X, ’ELEMENT  COORDINATES ’,/, 6X ,’ ELEMENT ’, 7X ,’ X-LEFT ’ , 
*7X, ’Y-LEFT’ ,6X, ’X-RIGHT ’ , 6X, ’Y-RIGHT’ ,//) 

DO  30  I»1,NELEMJ*2 

WRITE (6, 25)  I , ELEM( 1 , 1 ) , ELEM (1,2), ELEM (1,3), ELEM (1,4) 

25  FORMAT (9X, I  3, 5X, F 1 0 . 4, 3 ( 3X, FI  0 . 4 ) ) 

30  CONTINUE 
RETURN 
END 


81 


n  r> 
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c 

C  GSHATT 
C 

C  THIS  SUBROUTINE  MOVES  THE  CONNECTING  ELEMENT  TERMS  Of  THE  GS 
C  MATRIX  OVER  TO  THE  HS  MATRIX  AND  THEN  CALCULATES  GSHAT 
C 

SUBROUTINE  GSHATT 
implicit  real*8  (a-h,o-z) 

COMMON/GEOMET/NNODES , NELEMS , ORDER, EPS , NORMI , NORMJ , JAC, NODEC , 
*NODPIN 1 ,NODPIN2, NODE ( 5 ) , ELEM( 360,4 ) 

COMMON /VELOC IT/UVEL (240), WEL (240), XFORC (480), YFORC (480), 

*XOLD ( 480 ) , YOLD ( 480 ) 

COMMON/STRUCT/GE , XN  U , HS ( 1 92  0 , 1 92  0 ) , GS (1920,  1440), GS HAT (1920) 
INTEGER  ORDER 

REAL*8  NODE, NORMI, NORMJ, JAC 
C 

C  MOVE  THE  CONNECTING  ELEMENT  TERMS  OF  THE  GS  MATRIX  OVER  TO 
C  THE  HS  MATRIX 
C 

DO  10  1*1 , 8*NNODES 

DO  10  J* 1 , NNODES 

HS  ( 1 , 3*NNODES+J  )— GS  ( 1 , 2*NNODES+J  ) 

HS ( I , 7*NNODES+J )— GS (1,5 *NNODES+J ) 

10  CONTINUE 

CALCULATE  THE  GSHAT  MATRIX 
C  REMEMBER  THE  COLUMNS  OF  GS  CORRESPONDING  TO  THE  CONNECTING 
C  ELEMENTS  DO  NOT  COUNT 
C 

C  SUBTRACT  OLD  TRACTIONS  ( XOLD, YOLD )  FROM  NEW  TRACTIONS  ( XFORC, YFORC ) 
C  THIS  PRE-STRESSES  THE  STRUCTURE  TO  ACCOUNT  FROM  THE  DEFORMATION  AND 
C  STRESS  DUE  TO  THE  LAST  PASS 
C 

DO  20  1*1 , 8*N NODES 
GSHAT ( I )-0.0 
DO  20  J*1,2*NNODES 

GSHAT (I  )-GSHAT(I )+GS(I, J)*(X1  ^RC ( J )-XOLD( J ) )  + 

*  GS ( 1 , 3 ‘NNODES +J ) * ( YFORC ( J ) -YOLD ( J  )  ) 

20  CONTINUE 
RETURN 
END 
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C 
C 
C 
C 

C  THIS  SUBROUTINE  CALCULATES  THE  DIAGONAL  TERMS  O T  THE  H  MATRIX 
C 

c. . . 

SUBROUTINE  HDIAG 
implicit  real*8  (a-h,o-z) 

COMMON /GEOMET/NNODES , NELEMS , ORDER , EPS , NORM I , NORM J , JAC , NODEC , 
•NODP IN  1 , NODPIN 2 , NODE ( 5 ) , ELEM ( 36 0 , 4 ) 

COMMON/FLOW/H( 240 ,240),G(240,240),Q(240) , GHAT ( 240  ) 
integer  order 

real *8  J AC, NO RM I, NORM J, NODE 
C 

C  EXTERNAL  FLOW  (UNBOUNDED  INFINITE  DOMAIN) 

C  HII  —  SUM(HIJ)-M 
C 

DO  10  1“ 1 , NNODES 
SUM-0.0 

DO  3  J-1, NNODES 

IF(I.NE.J)  SUM  -  SUM+H{ I , J ) 

3  CONTINUE 

H ( I , I )— ( -SUM+1 . ) 

10  CONTINUE 

RETURN 
END 


83 


HGQOUT 


C 

C 

C 

C 

c 

c 

c 

c 


THIS  SUBROUTINE  OUTPUTS  THE  POTENTIAL  FLOW  MATRICES 
H,  G,  Q,  GHAT 


SUBROUTINE  HGQOUT 
implicit  real*8  {a-h,o-z) 

COMMON /GEOMET/NNODES , NELEMS , ORDER , EPS , NORMI , NORMJ , JAC , NODEC , 
•NODPifi  t .  NODPIN 2 , NODE  ( 5  )  ,  ELEM ( 360 , 4  ) 

COMMON/FLOW/H( 240 , 240),G(240,240),C(240) , GHAT ( 240 } 
real*8  NORMI,NO»1J,  JAC,  NODE 
INTEGER  ORDER 


12 

20 

21 

22 

31 

32 


WRITE (6, 1 ) 

FORMAT(//,22X, ' *****POTENTIAL  FLOW  MATRICES*****',/, 
22X, ' **  FOR  OUTER  SURFACE  BOUNDARY  **',//) 

WRITE (6, 2) 

FORMAT(5X,'H  MATRIX',/) 

DO  10  I  — 1 , NNCDES 

WRITE (6, 3 )  (H(I,J),J-1, NNODES  ) 

FORMAT ( /13(8(F10.4,5X) ,/) ) 

CONTINUE 
WRITE (6,11) 

FORMAT(//r5X, 'G  MATRIX',/) 

DO  20  1-1, NNODES 

WRITE{6, 12 )  (G(1,J),J-1, NNODES) 

FORMAT (/13(8(F10.4,SX) ,/) ) 

CONTINUE 
WRITE (6,21 ) 

FORMAT(//,5X, 'Q  MATRIX',/) 

WRITE (6,22)  (Q(I), 1-1, NNODES) 

FORMAT( 13(6 <F 10.4, 5X ) ,/) ) 

WRITE (6,31) 

FORMAT(//,5X, 'GHAT  MATRIX',/) 

WRITE (6,32)  (GHAT(I), 1-1, NNODES) 

FORMAT (13(8(F10.4,5X),/) ) 

RETURN 

END 
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C  HSDIAG 
C 

C  THIS  SUBROUTINE  CALCULATES  THE  DIAGONAL  TERM  OF  THE  HS  MATRIX 
C 

C  HI I  -  -SUM(HIJ) 

C 

C  FOR  THE  HXX  AND  HYY  TERMS  ONLY  BECAUSE  CXY-CYX-Q.  AND  THEREFORE 
C  THE  DIAGONAL  TERMS  OF  HXY  AND  HYX  ARE  EQUAL  TO  THEMSELVES 
C 

SUBROUTINE  HSDIAG 
implicit  real*8  (a-h,o-z) 

COMMON/GEOMET/N NODES , NELEMS .ORDER , EPS , NORM I , NORM J , J AC , NODEC , 
•NODPIN 1 , NOD PIN 2 , NODE (5),ELEM(360,4) 

COMMON/STRUCT/GE , XNU, HS ( 1 920,  192 0 ) ,GS (  1 92 0 ,  1  440  )  ,GSHAT(  192  0 ) 
real*8  NORMI ,  NORMJ,  JAC,  NODE 
INTEGER  ORDER 
C 

SUM  FOR  THE  OtTTER,  INNER  AND  CONNECTING  ROWS 

DO  10  1“ 1 , 3 ‘NNODES 
SUMXX -0.0 
SUMYY-0.0 
DO  5  J«1,3*NNODES 
IF(I.NE.J)  THEN 

SUMXX -SUMXX  +HS ( I , J ) 

SUMYY-SUMYY+HS ( 4 'NNODES  +  1 , 4  *NNODES+J ) 

END  IF 
5  CONTINUE 

HS ( I , I  )  =  -SUMXX 

HS ( 4 ‘UNODES+I , 4*NNODES+I !«*-SUMYY 
10  CONTINUE 
C 

C  CONTINUE  TO  SUM  FOR  THE  ADDITIONAL  SET  OF  CONNECTING  ROWS 
C 

DO  20  I-3*NNODES+ 1 , 4 *N NODES 
II-I-NNODES 
SUMXX -0 . 0 
SUMYY-0.0 
DO  15  J«1 , 3 *NNODES 
IF(II.NE.J)  THEN 

SUMXX-SUMXX+HS ( I , J ) 

SUMYY -SUMYY  +HS  (  4  *N NODES  + 1 , 4  *NNOD£S+ J ) 

END  IF 
15  CONTINUE 

HS(  I  ,  II  )=-SUMXX 

HS(4*NNODES+I , * *NNODES+I I )=-SUMYY 
20  CONTINUE 
RETURN 
END 
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C 

c 
c 
c 

C  THIS  SUBROUTINE  CALCULATES  THE  JACOBIAN  AND  OUTWARD  NORMAL 
C  FOR  THE  POINT  UNDER  CONSIDERATION 
C 

C  INPUT: 

C  J  -  ELEMENT  ON  WHICH  POINT  IS  LOCATED 

C  ETA  -  POINT  WHERE  JACOBIAN  AND  NORMAL  IS  CALCULATED 

C 

C  OUTPUT : 

C  JAC  -  JACOB  UN 
C  NORMI  -  X  COMPONENT  OF  NORMAL 
C  NORMJ  “  Y  COMPONENT  OF  NORMAL 
C 

SUBROUTINE  JACOBIAN ( J , ETA ) 
implicit  real*8  (a-h,o-r) 

COMMON/GEOMET/NNODES , NELEMS, ORDER, EPS , NORMI , NORMJ,  JAC , NoDEC , 
*NODPIN1,NODPIN2,NODE<5) ,ELEM{ 360,4) 
real *8  NORMI , NORMJ, JAC, NODE 
INTEGER  ORDER 
C 

C  ASSUME  FLAT  PANEL 
C 

DN1  —  .5 
DN3-.5 

DXDN«DN1*ELEM(J, 1 )+DN  3*ELEM( J , 3 ) 

DYDN»DN1*ELEM(J,2 )+DN3*ELEM( J,4 ) 

JAC«SyRT(DXDN*DXDN+DYDN *DYDN ) 

NORMI— DYDN/JAC 

NORMJ-DXDN/JAC 

RETURN 

END 
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Q* *****************************************************************  * 

c 

C  PRESSURE 

c 

C  THIS  SUBROUTINE  CALCULATES  THE  COEFFICIENT  OF  PRESSURE  (CP) 

C  AND  THE  EXTERNAL  FORCES  ( XFORC , YFORC )  ON  EACH  NODE 
C 

SUBROUTINE  PRESSURE 
implicit  real*8  (a-h,o-z) 

COMMON/GEOMET/NNODES , NELEMS , ORDER , EPS , NORM  I , NORMJ , JAC , NODEC , 
•NODPIN  1 , NODPIN2, NODE ( 5 ) , ELEM (360,4) 

COMMON/POTEN/VINF, ALPHA , VI , VJ , PINF , RHOINF , GAMMA, MINF2 , NODEP 
COMMON /VELOCIT/UVEL(  240  )  ,WEL(  240)  .XFORC  (480),  YFORC  (480)  , 
*XOLD(480 ) ,YOLD(480 ) 

REAL*8  JAC , NORM I , NORM  J , GAMMA , M INF  2 , NODE 
INTEGER  ORDER 

STORE  THE  OLD  VALUES  OF  TRACTIONS 

DO  10  I«1,2*NNODES 
XOLD ( I ) -XFORC (I ) 

YOLD( I )-YFORC ( I ) 

10  CONTINUE 

INITIALIZE  XFORC, YFORC 

DO  15  I  -  1 , 2*NNODES 
XFORC ( I ) —0 . 0 
YFORC ( I )-Q . 0 
15  CONTINUE 

CALCULATE  PRESSURE  AND  X,Y  COMPONENTS  OF  PRESSURE  FOR  OUTER  SURFACE 
XFORCE  -  NI*P,  YFORC  -  NJ*P 

DO  20  J-1, NELEMS 
DO  20  K-1, ORDER 

JJJ-J ‘ORDER- (ORDER-K) 

VEL-S£RT(UVEL( JJJ  )  *UVEL<  JJJ  )+WEL(  JJJ  )  *WEL(  JJJ  )  ) 

CP- 1 . - ( VEL/VINF ) • *2 
P-0. 5*RHOINF»VINF*VINF»CP 
CALL  JACOBIAN (J,NODE(K) ) 

XFORC ( J JJ ) -XFORC (JJJ)+P*( NORM I  ) 

YFORC (JJJ ) -YFORC ( JJJ )+P»( NORMJ  ) 

20  CONTINUE 

CALCULATE  PRESSURES  FOR  INSIDE  BOUNDARY  WHERE  INTERNAL  PRESSURE 
EQUALS  PRESSURE  AT  VENTED  NODE  ON  OUTER  SURFACE 

VELVENT-SQRT  (UVEL(  NODEP  )  *UVEL(  NODEP  )+WEL(  NODEP  )*WEL(  NODEP)  ) 
CPVENT- 1 . - ( VELVENT/VI NF ) *  *  2 
PIN«.5*RHOINF*VINF*VINF*CPVENT 
DO  30  J -NELEMS +1, 2 *NELEMS 
DO  30  K-1, ORDER 

JJJ -J ‘ORDER- ( ORDER-K ) 
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CALL  JACOBIAN ( J , NODE (K ) ) 

XFORC( JJJ )«XFORC< JJJ )+PIN*< -NOfWI ) 
Y  FORC <  J  J  J  )  -Y  FORC ( JJJ ) +P I N • ( -NORM J ) 
CONTINUE 
RETURN 
END 


PRESSURE  2/2 
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C 
C 
C 
C 

C  THIS  SUBROUTINE  CALCULATES  THE  BOUNDARY  CONDITION  VECTOR, Q 
C  USING  Q  «  -n  .  Vinf 
C 

C  NOTE  THAT  SINCE  THE  OUTWARD  NORMALS  TO  THE  FLUID  FLOW  POINT  INTO  THE 
C  BODY  THE  NORMAL  CALCULATED  IN  JACOBIAN  SUBROUTINE  ALREADY  HAS  THE 
C  PROPER  SIGN  (-N) 

C 

SUBROUTINE  QFLOW 
implicit  real*8  (a-h,o-z) 

COMMON/GEOMET/NNODES , NELEMS , ORDER, EPS , NORM  I , NORM J , JAC , NODEC , 
•NODPIN 1 , NODPIN  2 , NODE ( 5 ) , ELEM (360,4) 

COMMON/POTEN/VINF , ALPHA , VI , VJ , PINF , RHOINF , GAMMA , MINF2 , NODEP 
COMMON/FLOW/H ( 240, 240 ) ,G(240,240 ) ,Q(240 ) , GHAT (240 ) 

INTEGER  ORDER 

real *8  JAC, NORM I, NORM J, GAMMA, MINF 2, NODE 
C 

C  EXTERNAL  FLOW  PROBLEM 
C  Q “NORMAL* VINF 
C 

DO  5  I“1 , NNODES 
Q( I )“0. 0 

5  CONTINUE 

DO  10  J-1, NELEMS 
DO  10  L-1, ORDER 

JJ-J*ORDER-(ORDER-L) 

CALL  JACOBIAN ( J,NODE(L)  ) 

Q(JJ  J“Q(JJ  )+(NORMI*VI+NORMJ*VJ  ) 

10  CONTINUE 

RETURN 
END 
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C 
C 
C 
C 

C  THIS  SUBROUTINE  CALCULATES  THE  H  AND  G  MATRICES  FOR  THE 
C  ELASTIC  STRUCTURE 
C 

SUBROUTINE  SFORM 

implicit  real'8  (a-h.o-t) 

COMMON /GEOMET/NNODES , NELEMS , ORDER , EPS , NORMI , NORMJ , JAC , NODEC , 
'NODPIN 1 , NODPIN2 , NODE ( 5 ) , ELEM (360,4) 

COMMON/STRUCT/GE , XN U , HS ( 1 92  0 ,  1 92  0  ) , GS (  1 92  0 ,  1 44 0  ) , GSHAT <  1 92  0 ) 

COMMON /POSIT ION /XI ( 1080) ,YI< 1080 ) 

COMMON/VELOCIT/UVEH  240  )  ,  WEL(  240  )  ,XFORC  ( 480 ) , YFORC ( 480 ) , 

•XOLD (480), YOLD (480) 

REAL*8  ETA( 12) ,WT( 12) ,ETALN< 10 ) , WTLN ( 1 0  )  , ETAS ( 8  )  ,WTS(8 ) , 

*ETA2( 10) ,WT2( 10) 

REAL'S  SHAP( 5 ) , SHAP2 ( 5 ) , SHAPLN ( 5 ) , SHAPLN2 ( 5 ) ,SHAP3 ( 5 ) , 

•  SHAP4 ( S ) 

REAL*8  NODE . JAC , JAC 1 , JAC 2 , JAC  3 , JAC 4 , JACLN , JACLN 2 , NORM I , NORMJ , 
'NORMI 1 , NORMJ 1 . NORTI 2, NORMJ 2 

INTEGER  ORDER 
C 

C  ETA  AND  WT  ARE  12TH  ORDER  GAUSSIAN  QUADRATRUE  VALUES  AND  WEIGHTS 
C 

DATA  ETA  (  1  )  ,  ETA  (  2  )  ,  ETA  (  3  )  ,  ETA  {  4  )  ,  ETA  (  5  )  ,  ETA  (  6  )  ,  ETA  (  7  )  , 

'ETA ( 8 ) , ETA ( 9 ) , ETA ( 10 ) , ETA (  1 1 ) , ETA (  12 ) /-. 98 1 56063d0, -. 904 1  1 726d0 , 
*-. 76990267d0,-. 5873 179Sd0, -.3678314980,-. 12523341d0, . 12523341d0, 
*. 367831 49d0, .5873 1795d0, .  76990267d0 , . 904 11 726d0, .981 56063d0/ 

DATA  WT( 1 ) , WT( 2 ) ,WT(3 ) , WT(4 ) , WT( 5 ) , WT( 6 ) , WT( 7 ) , WT(8 ) , WT( 9 ) , 

•WT (  10 ) ,WT ( 1 1 ),WT(  12 )/.0471753dO, .  10693933d0, .  16007833d0, 

• . 203 16743d0 , .23349254d0,  .24914705d0, . 249 1 4705d0 ,. 233492 54d0 , 

* . 203 16743d0, . 16007833d0, . 10693933d0, . 04717534d0/ 

C 

C  ETALN  AND  WTLN  ARE  10TH  ORDER  LOGARITHMIC  GAUSSIAN  QUADRATURE 
C  VALUES  AND  WEIGHTS 
C 

DATA  ETALN ( 1 ) , ETALN ( 2 ) , ETALN ( 3 ) , ETALN (4  )  .ETALN ( 5  )  .ETALN (6  )  , 
'ETALN ( 7), ETALN (8),  ETALN (9),  ETALN (  10 )/. 00904259d0, .05397105d0, 

*. 13531 134d0, .24705169d0, . 3802 1 1 7 1 dO , . 52 3791 59d0 , . 66577472d0 , 

'. 7941 90 19d0, . 89816 102d0, .96884798d0/ 

DATA  WTLN ( 1 ) , WTLN (2 ) , WTLN  < 3 ) , WTLN ( 4 ) , WTLN ( 5 ) , WTLN (6 ) , WTLN ( 7 ) , 
'WTLN (8), WTLN (9),  WTLN (  10  )/ .  1 20 95474d0 , .  1 86363 1 0dO 1 9566066d0 , 

*.  1 735772 3d0, .  13569597d0,  .093  64  70  8d0, . 0 5578793 8d0 , .02715989d0, 

*. 00951 5 1992d0, .00 16381 586d0/ 

C 

C  ETAS  AND  WTS  ARE  8TH  ORDER  GAUSSIAN  VALUES  AND  WEIGHTS  FOR 
C  1/R  SINGULAR  INTEGRALS 
C 

DATA  ETAS ( 1 ) , ETAS (2 ) .ETAS (3 ) .ETAS (4 ) .ETAS (5 ) , ETAS ( 6 ) , ETAS ( 7 ) , 
*ETAS(8 )/-. 00324250d0, . 053490 77d0, . 17782733d0, .350?1788d0, 

54581 95sd0, . 73342708d0 , . 88498305d0 , .9774543880/ 
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DATA  WTS (  1  ) , WTS ( 2  ) ,  WTS  (  3  )  ,  WTS ( 4 ) , WTS ( 5 ) , WTS ( 6 ) , WTS ( 7  ) , WTS ( 8  ) 

*/-3 . 93063681 dO, 1 . 7374063 1 dO , . 8S7634S4d0 , . 53836077d0 , 

*. 35975 759d0, . 23726769d0, . 141462 28 dO, . 05874 761d0/ 

C 

C  ETA 2  AND  WT2  ARE  10TH  ORDER  GAUSSIAN  QUADRATURES 
C 

DATA  ETA2 ( 1 ) , ETA 2(2), ETA2 ( 3 ) , ETA 2 ( 4 ) , ETA 2 ( 5 ) , ETA 2 ( 6 ) , ETA 2 ( 7  ) , 
*ETA2 (8 ) , ETA2 ( 9 ) , ETA2 (10) /-. 97390653 dO,-. 8650 63 37d0,-. 67 94 0957 dO, 
*-.43339539d0,-. 14887434d0, . 14887434d0, . 43339539d0 , . 67940957d0 , 
*.86506337d0, .97390653d0/ 

DATA  WT2 ( 1 ) , WT2 ( 2 ) ,WT2 ( 3 ) , WT2 ( 4 ) , WT2 ( 5 ) , WT2 (6 ) ,  WT2  ( 7 ) , WT2 ( 8 ) , 
*WT2 ( 9 ) , WT2 (10)/. 0666 7134 dO, . 149451 3 5d0, . 2 1 908636d0, . 269266 72 dO , 
*.295S2423d0, .29552423d0,  .26926672d0, .21908636d0, . 149451 35d0, 

*. 06667 134d0/ 

DATA  PI/3. 14159265d0/ 

C 

C  INITIALIZE  MATRICES 
C 

DO  15  I-1,8*NNODES 
DO  10  J-1,6*NNODES 
GS(I,J)-0.0 
10  CONTINUE 

DO  15  J-1,8*NNODES 
HS(I,J)-0.0 
15  CONTINUE 
C 

C  BEGIN  LOOPS  ON  NODE  I  AND  ELEMENT  J 
C 

DO  1000  111-1,4 

DO  1000  II - 1 , HELENS 
DO  1000  K-1, ORDER 

I-( III -1  )  *NNODES+I I ‘ORDER- (ORDER-K ) 

IFUII.EQ.4)  THEN 

INODE-2 *NNODES+II*ORDER- (ORDER-K) 

ELSE 

INODE -I 
END  IF 

DO  1000  L-1,4 

IF( III . EQ. 4 )  THEN 
J-NELEMS‘(L-1 )+II+1 

IF(L.EQ.4. AND.II.NE.NELEMS )  J-2‘NELEMS+II 
IF(II.EQ.NELEMS )  J -NELEMS * (L-1  )  + 1 
IF( 1 1 . EQ. NELEMS . AND • L. EQ. 4 )  J-3 ‘NELEMS 
ELSE 

J-NELEMS* ( L-1 )+II 

IF  (  L.  EQ.  4  )  J-2 ‘NELEMS'*’ ( II  - 1  ) 

IF(L.EQ.4.AND.II.EQ. 1)  J-3*NELEMS 
END  IF 
C 

C  DETERMINE  IF  THE  INTEGRATION  SHOULD  BE  SINGULAR 
C 

C  I.E.  IS  THE  NODE  I  SOMEWHERE  ON  THE  ELEMENT  J 
C  ISING  -  0  MEANS  THE  NODE  I  IS  NOT  ON  ELEMENT  J 
C  ISING  -  1  MEANS  THE  NODE  I  IS  THE  FIRST  NODE  ON  ELEMENT  J 
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C  ISING  -  2  MEANS  THE  NODE  I  IS  THE  SECOND  NODE  ON  ELEMENT  J 
C  ISING  “  3  ETC 
C 

ISING-0 

DO  20  LL-1, ORDER 

NODES -J ‘ORDER- ( ORDER-LL ) 

IF ( INODE. EQ. NODES )  ISING-LL 
20  CONTINUE 

IF( ISING. EQ. 0 )  THEN 

c  «•««««  t  i  «  t  ••»«*«««  t  >.  t  t  »•  « 

C  THE  INTEGRATION  IS  NOT  SINGULAR 
C  PERFORM  1 2TH  ORDER  GAUSSIAN  QUADRATURE 

£••*••••••••••••••••»»•••••*••• 

DO  50  *3C«1 ,  12 
C 

C  CALCULATE  THE  X  AND  Y  COORDINATES  OF  THE  QUADRATURE  POINT 
C  AND  CALCULATE  THE  DISTANCE  FROM  NODE  I  TO  QUAD  POINT 
C 

ETA  A  —  1 . 

ETAC-1. 

ETAQ-ETA ( KK ) 

SHAP ( 1 )» ( ETAQ-ETAC ) / ( ETAA-ETAC ) 

S  HAP ( 2 ) - ( ETAQ -ETAA  J / ( ETAC-ETAA  ) 

X«SHAP< 1 ) *EL£M( J , 1 ) +SHAP ( 2 ) *ELEM ( J , 3 ) 
Y-SHAP1 1 ) *ELEM ( J , 2  )  +SHAP ( 2 ) *ELEM( J , 4 ) 
RI-X-XI( INODE) 

RJ-Y-YI  ( INODE) 

R-SQRT( RI »RI+RJ*RJ ) 

CALCULATE  U*  AND  P* 

C 

CALL  JACOBIAN (J, ETA ( KX) ) 

IF(L.EQ. 1.0R.L.EQ.3)  THEN 
NORM  I— NORM  I 
NORMJ--NORMJ 
END  IF 

TERM 1*( -1 . )/(8. *PI *GE • ( 1 . -XNU ) ) 

C 

C  NOTE  THAT  ON  THE  LEFT  CONNECTING  ELEMENT  (L-4)  Qij— Qji 
C  SO  THAT  THE  GS  TERM  SHOULD  BE  NEGATIVE 
C 

IF (L. EQ. 4 )  TERM  1 “-TERM  1 
USTXX-TERM1*<  < 3 . -4 . *XNU ) *dlog( 1./R) 

*  ♦( RI/R ) “2  ) 

USTXY-TERM1*((RI*RJ)/(R»R) ) 

USTYXHISTXY 

USTYY-TERM1*( ( 3 . -4 . »XNU ) ‘dlog ( 1./R) 

*  +  ( RJ/R)  “2  ) 

DRDN  -  ( NORM  I  *  RI -t-NORM  J  *R  J  )  /  R 
TERM2-(-1. )/(4.*PI*( 1 • -XNU ) *R ) 
PSTXX-TERM2»(DRDN*l ( 1 . -2. ®XNU)+ 

*  2. * ( RI/R) **2 ) ) 
PSTXY«TERM2*(DRDN*2.»(RI»RJ)/(R*R) 

*  -( 1.-2. *XNU ) * ( (  RI  ‘NORMJ-RJ ‘NORM I )/R) ) 
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PSTYX*TERM2*(  DRDH *2 . * ( R1 *RJ )/( R‘R ) 

*  -( 1.-2.*XNU)*( (RJ*NORMI-RI*NORMJ)/R) ) 
PSTYY-TERM2*(DRDN*( ( 1.-2.*XNU)‘ 

*  2.*(RJ/R)“2  )  ) 

C 

C  CALCULATE  THE  TERMS  OF  H  AND  G  MATRICES 
C  ROW  I  CORRESPONDS  TO  THE  X  DEFORMATIONS  OF  NODE  I 
C  ROW  NNODES+I  CORRESPONDS  TO  THE  Y  DEFORMATIONS  OF  NODE  I 
C  COLUMN  JJ  CORRESPONDS  TO  THE  X  TERMS 
C  COLUMN  NNODES  +  JJ  CORRESPONDS  TO  THE  Y  TERMS 
C 

C  FOR  EXAMPLE, 

C 

C  -  - 

C  f  PXX*1  PXX*2  .  .  PXX ‘ORDER  PXY*1  PXY*2  .  •  PXY ‘ORDER  1 

C  *  PYX* 1  PYX*2  .  .  PYX‘ORDER  PYY*1  PYY*2  .  .  PYY ‘ORDER  * 

C 

CALL  SHAPE (ORDER, ETA(KK  )  ,SHAP) 

DO  35  LL-1, ORDER 

JJ-J ‘ORDER- (ORDER-LL ) 

HS(I,JJ)-HS(I,JJ)‘PSTXX*SHAP(LL)*JAC*WT(KK) 
HS ( 1 , 4 ‘NNODES+JJ  )  -HS (1,4 *NNODES‘JJ ) 

*  +PSTXY  *SHAP ( LL  )  *JAC‘WT( KX ) 

HS ( 4‘NNODES+I , JJ ) -HS ( 4*NNODES+I , JJ ) 

*  +PSTYX*SHAP(LL)*JAC*WT(KK) 

HS <  4‘NNOOES+I , 4 ‘NNODES +J J )- 

*  HS ( 4 ‘NNODES+1 , 4‘NNODES+JJ ) 

*  +PSTYY*SHAP(LL)*JAC*WT(KX) 

GS ( I , JJ )-GS ( I , JJ )‘USTXX‘SHAP(LL ) *JAC*WT ( KK ) 
GS  ( I  ,  3 ‘NNODES ■‘JJ  )-GS  ( 1 , 3‘NNODES+JJ  ) 

*  +USTXY  *SHAP( LL  )  *JAC*WT ( KX ) 

GS ( 4 ‘NNODES +1 , JJ ) -GS ( 4 ‘NNODES ‘I , JJ ) 

*  +USTYX*SHAP (LL )  ‘JAC*WT(KX ) 

GS ( 4 ‘NNODES  +  1 , 3  *NNODES+JJ )- 

*  GS ( 4 ‘NNODES+1 , 3*NNODES+JJ ) 

*  +USTYY*SHAP(LL) *JAC‘WT(KX) 

35  CONTINUE 

50  CONTINUE 

ELSE 

(:•••••••••••••• 

C  INTEGRAL  IS  SINGULAR  INSIDE  THE  ELEMENT 
C 

C  CALCUALTE  HS  AND  GS  MATRICES  SEPARATELY  BECAUSE 
C  GS  IS  LOG ( 1/R )  SINGULAR  AND  HS  IS  (1/R)  SINGULAR 
C 

C  CALCULATE  GS  FIRST 
C 

C‘‘‘‘‘‘*‘‘‘‘‘‘‘‘‘*‘‘‘‘‘‘‘‘‘‘‘‘‘‘ 

SINGUL-NODE ( I SING ) 

X2X 1-ELEM( J , 3 )-ELEM ( J , 1  ) 

Y2Y 1 -ELEM( J , 4 ) -ELEM ( J, 2 ) 
PLENO-SORT<X2X1*X2X1+Y2Y1*Y2Y1 ) 

DO  350  XX- 1,10 
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325 

350 


CALL  SHAPE (ORDER, ETA 2 (KK ) ,  SHAP ) 

CALL  JACOBIAN (J, ETA2 (KK  ) ) 

JAC2-JAC 

ET«S INGUL- ( S INGUL* 1 . ) *ETALN ( KK ) 

CALL  SHAPE (ORDER, ET, SHAPLN  ) 

CALL  JACOB IAN ( J , ET ) 

JACLN-JAC 

ET“S INGUL* ( 1 . -SINGUL ) *ETALN (KK ) 

CALL  SHAPE (ORDER, ET, SHAPLN2 ) 

CALL  JACOB IAN ( J , ET ) 

JACLN2-JAC 

ET-SINGUL-( 1 .+SINGUL)*.5*(ETA2(KK)+1 . ) 

CALL  SHAPE! ORDER, ET,SHAP3) 

CALL  JACOB IAN (J,ET) 

JAC3-JAC 

ET-SINGUL*( 1. -SINGUL )*.5*(ETA2( KK)  +  1  .  ) 

CALL  SHAPE (ORDER, ET,SHAP4) 

CALL  JACOB IAN (J,ET) 

JAC4-JAC 

DO  325  LL-1, ORDER 

JJ “J ‘ORDER- (ORDER-LL ) 

TERM«(-1. )/(  8 . *PI *GE* ( 1 . -XNU ) ) 

IF(L.EQ.4  )  TERM— TERM 
TERMX*(X2X 1/PLENG ) **2*TERM 
TERMY-(Y2Y1/PLENG)**2»TERM 

TERM3-(3.-4.*XNU)*dlog(ABS(PLENG*(  1.+SINGUD/2.  )  ) 
♦TERM 

TERM4*( J . -4 . *XNU ) *  tiiog ( ABS i PLENG* ( SINGUL- 1 . )/2. ) ) 
•TERM 

TERM  1 - ( 3 . -4 . *XNU ) ‘TERM 
TERM2»( 3. -4. *XNU) ‘TERM 
GS(I,JJ)-GS(I,JJ)+ 

SHAP ( LL ) »TE RMX * JAC2 *WT2 ( KK  ) - 

SHAP3'LL)*T£RM3*JAC3»(  1  -  *SINGUL ) * . 5*WT2 (KK ) - 
SHAP4 ( LL ) *TERM4  *JAC4* ( 1 . -SINGUL )*. S«WT2 ( KK ) + 
SHAPLN ( LL ) ‘TERM  1  * J ACLN • (1 . +S INGUL  )  *WTLN  ( KK )  + 
SHAPLN 2 ( LL) ‘TERM 2  *JACLN2* ( 1 . -SINGUL ) *WTLN ( KK ) 

GS (1,3  *NNODES  + J  J )-GS(I , 3*NNODES+JJ )  + 

SHAP ( LL ) *X  2X 1 *Y  2Y 1  * J AC  2  *WT2 ( KK ) *TE RM/PLENG  *  *2 
GS(4*NNODES*I,JJ)-GS(4*NNODES+I ,JJ)* 

SHAP( LL ) *X2X 1 *Y  2Y 1  *JAC2*WT2  (  KK )  *TE RM/PLENG* *2 
GS  <  4*NNODES*I , 3*NNODES*JJ  )- 
GS(4*NNODES+I,3*NNODES+JJ  )  + 

SRAP( LL ) ‘TERM Y*JAC2*WT2 ( KK ) - 

SHAP 3 ( LL ) ‘TERM 3  *JAC3* ( 1 . *SINGUL ) * . 5 *WT2 ( KK ) - 
SHAP4 ( LL ) *TERM4*JAC4* ( 1 . -S INGUL  )*. 5*WT2(KK )  + 
SHAPLN ( LL } ‘TERM  1 *JACLN  * ( 1 . +SINGUL ) *WTLN (KK)* 
SHAPLN  2 ( LL ) *TERM2  * JACLN  2* ( 1 . -SINGUL ) *WTLN (KK ) 
CONTINUE 
CONTINUE 


C  CALCULATE  THE  HS  MATRIX 
C  FOR  SINGULAR  PANELS  PXX-PYY-0. 
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DO  450  KX-1 , 6 

TERM«( 1 . -2.*XNU)/(4. ‘Pl*( 1 . -XNU) ) 

BT-(-1. -SINGUL) • ETAS ( KK ) +S INC UL 
CALL  JACOBIAN (J,ET) 

J  AC  1  “J  AC 
NORM  1 1  -NORM1 
NORMJ1-NORMJ 

IF(L.EQ. 1.0R.L.EQ.3)  THEN 
NOIMI 1 — NORMI 1 
NORM J 1 “-NORM J  1 
END  IF 

PXY 1 “TERM* ( (X2X1/PLENG)‘NORMJ1-(Y2Y1/PLENG)*NORMI1 ) 
PYX1 “TERM* ( ( Y  2Y 1 /PLENG ) *NORM1 1- (X  2X 1 /PLENG ) ‘NORMJ 1 ) 
CALL  SHAPE (ORDER, ET.SHAP) 

ET-(  1. -SINGUL) ‘ETAS (KK)+SINGUL 
CALL  JACOB I AN ( J , ET ) 

JAC2-JAC 

NORMI2-NORM1 

NORMJ2-NORMJ 

IF(L.EQ.  1.0R.L.EQ.3)  THEN 
NORMI2— NORMI2 
NORMJ2— NORMJ2 
END  IF 

PXY2«TERM‘< < X2X 1 /PLENG )*NORMJ2-(Y2Yl /PLENG) ‘NORM  12 ) 
PYX2«TERM*(  ( Y  2Y  1 /PLENG  )  *NORMI  2- ( X2X  1  /PLENG  )  *NOFMJ2  ) 
CALL  SHAPE (ORDER, ET,SHAP2) 

DO  425  LL-1, ORDER 

JJ-J ‘ORDER- ( ORDER-LL ) 

HS ( 1 , 4 *NNO DES ♦ J J ) -HS ( 1 , 4 ‘NNODES+JJ  )  - 
‘  SHAP(LL)*PXY1*JAC1‘(2./PLENG)‘WTS(KK)+ 

*  SHAP2 (LL ) *PXY2 ‘JAC2* ( 2 . /PLENG ) *WTS ( KX ) 
HS(4‘NNODES‘I , JJ )-HS ( 4 ‘NNODES+I , JJ  )- 

*  SHAP ( LL ) ‘PYX 1  * JAC 1 • ( 2 . /PLENG ) ‘WTS (KK 1  + 

*  SHAP2 ( LL ) *PYX2‘JAC2* ( 2 . /PLENG ) *WTS ( KK ) 

425  CONTINUE 

450  CONTINUE 

C 

C  CALCULATE  F( SINGULARITY )  LOG  TERMS  AND  ADD  TO  THE  HS  TERMS 
C 

TERM»( 1.-2. *XNU)/(4. *PI‘( 1 . -XNU)  ) 

CALL  SHAPE (ORDER, SINGUL, SHAP3  ) 

CALL  JACOBIAN(J, SINGUL) 

IF(L.EQ. 1.0R.L.EQ.3 )  THEN 
NORM I --NORM I 
NORMJ--NORMJ 
END  IF 

PXY -TERM* ( ( X2X 1/PLENG)*NORMJ-(Y2Y 1 /PLENG  ) ‘NORM I ) 
PYX-TERM* ( (Y 2Y1 /PLENG ) ‘NORMI- (X2X1 /PLENG) ‘NORMJ ) 

DO  475  LL-1, ORDER 

JJ-J ‘ORDER- (CRDER-LL ) 

HS(I , 4 ‘NNODES+JJ )»HS ( I , 4*NNODES+JJ  )- 

*  SHAP3 ( LL ) *PXY  *JAC* ( 2 . /PLENG ) *dlog( ABS ( - 1 . -SINGUL ) )  + 

*  SHAP3 (LL ) ‘PXY ‘JAC* ( 2 . /PLENG ) *dlog<  ABS ( 1 . -SINGUL ) ) 

HS ( 4 ‘NNODES+I ,JJ ) -HS ( 4 ‘NNODES+I , JJ )- 
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SFORM  7/7 

SHAP3  ( LL  )  *PYX*JAC* (  2  .  /PLENG  )  «dlog<  ABS( -  1  .  -S  INGUI. )  )  + 
SHAP3 ( LL ) *PYX  *JAC* ( 2 . /PLENG ) *dlog( ABS ( 1.-SINGUL) ) 
CONTINUE 
END  IF 

CONTINUE 

RETURN 

END 
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c 

C  SHAPE 
C 

C  THIS  SUBROUTINE  CALCULATES  THE  SHAPE  FUNCTIONS 
C  USING  LAGRANG IAN  POLYNOMIAL  SHAPE  FUNCTIONS 
C 

C  INPUT: 

C  SORDR  -  ORDER  OF  SHAPE  FUNCTION 

C  TERM  -  LOCATION  WHERE  SHAPE  FUNCTION  IS  TO  BE  CALCULATED 
C 

C  OUTPUT : 

C  SH  »  SHAPE  FUNCTION  ARRAY 
C 

C*.* 

SUBROUTINE  SHAPE ( SORDR, TERM, SH ) 
implicit  real'8  (a-h,o-z) 

COMMON /GEOM ET/N NO DE S , NELEMS ,ORDE  R , EPS , NORM  I , NORM J , JAC , NODEC , 
•NODPIN 1 , NODPIN2 , NODE ( S ) ,ELEM(360,4) 
real's  SH(S  )  , NORM I , NORM J, JAC, NUMER, NODE 
INTEGER  ORDER, SORDR 

FOR  CONSTANT  FIT 

IF (SORDR. EQ. 1 )  THEN 
SH ( 1  )«1. 

ELSE 

CALCULATE  HIGHER  ORDER  SHAPE  FUNCTIONS 
C 

DC.  60  1-1,  SORDR 
DENOM-  1 . 

NUMER-1 . 

DO  55  K-1, SORDR 
IF(I.NE.K)  THEN 

DENOM«( NODE ( I )-NODE(K) ) *DENOM 
NUMER«(TERM-NODE (K) ) *NUMER 
END  IF 
55  CONTINUE 

SH ( I ) “NUMER /DENOM 
60  CONTINUE 

END  IF 
RETURN 
END 


SINFO 


C 

c 
c 
c 

C  THIS  SUBROUTINE  READS  IN  THE  INFORMATION  NEEDED  TO 
C  CALCULATE  THE  STRUCTURE  MATRICES 
C 

C  INPUT: 

C  GE  -  SHEAR  MODULUS 
C  XNU  -  POISSON  MODULUS 
C 

c..******»»....*.. 

SUBROUTINE  SINFO 
implicit  real*8  (a-h,o-z) 

COMMON /STRUCT/GE, XNU, HS( 1 92 0 ,  192 0  !  ,GS ( 192 0 ,  1 440 ) , GSHAT ( 192  0 ) 
C 

c 

READ ( 5 ,  10  )  GE , XNU 
10  FORMAT (2F 10. 4 ) 

RETURN 

END 
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SOLVE 


c 
c 

C  THIS  SUBROUTINE  SOLVES  THE  LINEAR  SYSTEM  OF  EQUATIONS 
C  BY  THE  GAUSS  ELIMINATION  METHOD  PROVIDING  FOR 
C  INTERCHANGING  ROWS  WHEN  ENCOUNTERING  A  ZERO  DIAGONAL 
C  COEFFICIENT 
C 

C  A:  SYSTEM  MATRIX 

C  B:  ORIGINALLY  THE  INDEPENDENT  COEFFICIENTS  ( RHS ) 

C  AFTER  SOLUTION  IT  CONTAINS  THE  VALUES  OF 

C  SYSTEM  UNKNOWNS 

C  Ns  ACTUAL  NUMBER  OF  UNKNOWNS 
C  D:  DETERMINANT 
C 

c**»* ******* . . 

SUBROUTINE  SOLVE ( A,B,D,N,NDIM) 
implicit  real*8  (a-h,o-z) 

REAL*8  A ( NDIM, NDIM) , B (NDIM) 

N1-N-1 

DO  10U  K-1.N1 
K1-K+1 
C-A(K,K) 

IF( ABS (C ) -0 . 00000 1 dO )  1,1,3 
1  DO  7  J-K1,N 

TRY  TO  INTERCHANGE  ROWS  TO  GET  NON  ZERO  DIAGAONAL 

IF( ABS ( A( J ,K)) -0.000001  HO)  7,7,5 

5  DO  6  L-K,N 

C»A(K,L) 

A(K,L)-A(J,L) 

6  A(J,L)-C 
C-B(K) 

B(K)-B(J) 

B(J)-C 

C-A(K,K) 

GO  TO  3 


7  CONTINUE 

8  WRITE (6, 2)  K 

2  PORMAT< ' ••••♦SINGULARITY  IN  ROW',15) 

D-0. 


GO  TO  300 
C 

C  DIVIDE  ROW  BY  DIAGONAL  COEFFICIENT 
C 

3  C=A(K,K) 

DO  4  J«K1,N 

4  A(K,J)-A{K, J)/C 
B(K)-B(K)/C 

C 

C  ELIMINATE  UNKNOWN  X(K)  FROM  ROW  I 
C 
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SOLVE  2/2 


DO  10  I-K1,N 
C-A(I,K) 

DO  9  J-K1,N 

9  A(I,J)-A(I,J)-C*A(K,J) 

10  B(I )-B(  I }-C*B(K) 

100  CONTINUE 

COMPUTE  LAST  UNKNOWN 

IF  ( ABS ( A ( N , N) )-0. 000001 dO )  101,101,102 

101  WRITE (6,2)  N 
D-0.0 

GO  TO  300 

102  B(N)«B(N)/A(N,N) 

APPLY  BACKSUBSTITUT ION  TO  COMPUTE  REMAINING  UNKNOWNS 

103  DO  200  L-1.N1 

K-N-L 

K1-K+1 

DO  200  J-K1.N 
200  B(K)»B(K)-A(K,J)*B(J) 

COMPUTE  VALUE  OF  DETERMINANT 


D-1 . 

DO  250  1-1, N 
250  D-D*A( 1,1) 
300  RETURN 
END 
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SOUT 


c 
c 
c 
c 

C  THIS  SUBROUTINE  PRINTS  OUT  THE  ELASTIC  STRUCTURE  MATRICES, 

C  HS,  GS  AND  GSHAT 
C 

SUBROUTINE  SOUT 
implicit  real*8  (a-h,o-z) 

COMMON/GEOMET/NNODES , NELEMS , ORDE  R , EPS , NORM I , NORM J, JAC , NODEC , 
•NODPIN 1 , NODPIN  2 , NODE ( 5 ) , ELEM (360,4) 

COMMON/STRUCT/GE , XNU, HS (  1920,  192 0  ) , GS (  192 0 ,  1440),GSHAT<  1920 ) 
real*8  NORMI , NORMJ , J AC, NODE 
INTEGER  ORDER 


WRITE ( 6 , 1 ) 

1  FORMAT ( // , 22X, ' *  **  “STRUCTURAL  MATRIX  INFORMATION*****’,//) 
WRITE (6,2) 

2  FORMAT ( 5X, 'HS  XX  MATRIX',/) 

DO  10  I-1,NNODES‘4 

WRITE  (6, 3)  (HS(I,J),J-1,NNCDES*3) 

3  FORMAT(35( 13(F6.4, 3X) ,/) ) 

10  CONTINUE 
WRITE (6,11) 

11  FORMAT ( 5X, 'HS  XT  MATRIX',/) 

DO  13  I  - 1 , NNODES *4 

WRITEI6, 12)  (HS ( I , 4*NNODES+J ) ,  J-1  , NNODES *3 ) 

12  FORMAT (35(13(F6.4,3X),/)) 

13  CONTINUE 
WRITE (6, 120 ) 

120  FORMAT (5X, 'HS  TX  MATRIX',/) 

DO  150  I« 1 , NNODES *4 

WRITE ( 6 , 115)  (HS( 4*NNODES+I , J ) , J-1 , NNODES *3 ) 

115  FORMAT (35(13(F6.4,3X) ,/)) 

150  CONTINUE 

WRITE (6,220  ) 

220  FORMAT ( 5X , ' HS  YT  MATRIX',/) 

DO  250  1“ 1 , NNODES *4 

WRITE (6, 21 5 )  ( HS ( 4  *NNODES  +  I , 4 ‘NNODES +J ) , J-1 , NNODES *3 ) 

215  FORMAT (35( 13(F6.4, 3X) ,/) ) 

250  CONTINUE 

WRITE (6, 14  ) 

14  FORMAT( //, 5X , 'GS  MATRIX',/) 

DO  20  1“ 1 , NNODES  *8 

WRITE ( 6 , 15 )  (GS( I , J ) , J-1 .NNODES *6 ) 

15  FORMAT ( 70 ( 1 3 ( F6 >  4 , 3X) ,/) ) 

20  CONTINUE 

WRITE (6,22  ) 

22  FORMAT!//, 5X, ’GSHAT  MATRIX* ,/) 

WRITE ( 6,23)  (GSHAT (I), I- 1,8 ‘NNODES) 

23  FORMAT! 95 (13<F7.4,3X),/)) 

RETURN 

END 
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SO  ITT  PUT 


C 
C 
C 
C 

C  THIS  SUBROUTINE  OUTPUTS  THE  DEFORMATION  MATRIX,  U 
C  NOTE  THAT  AFTER  THE  SOLVE,  THE  VECTOR  GSHAT  NOW  CONTAINS  THE 
C  SOLUTION  DEFORMATIONS 
C 

C*** ••••••••••••••••••••••••••••»••. 

SUBROUTINE  SOUTPUT ( ITER ) 
implicit  real*8  (a-h,o-i) 

COMMON /GEOMET/NNODES , NELEMS , ORDER, EPS , NORMI , NORM J , JAC , NODEC , 
•NODPIN  1 , NOD PIN 2, NODE (5),ELEM(360,4) 

COMMON/STRUCT/GE, XNU, HS ( 1920, 1920 ) ,GS ( 1920, 1440) ,GSHAT( 1920 ) 
real*8  NORMI, NORMJ, JAC, NODE 
INTEGER  ORDER 


WRITE ( 6 ,  1  )  GE , XNU, NODPIN 1 , NODPIN 2 , ITER-1 

1  FORMAT( //, 22X, ' *****DEFORMATION  MATRIX** ***' ,/, 

*1  OX, 'SHEAR  MODULUS  «  ', FI 0 . 4 , 5X , 'POISSON  MODULUS  »  ',F10.4,/, 
•10X, 'PINNED  NODES  ARE  -  ', 21 3/ , 1 0X ,' ITERATION  -  *,12./) 

WRITE (6, 2) 

2  FORMAT<//,SX, 'U  MATRIX* ,/) 

DO  10  I-1,4*NNODES 

WRITE (6,3)  GSHAT ( I ) ,GSHAT(4*NNODES+I  ) 

3  FORMAT ( F20. 8, 2X, F2  0 . 8 ) 

10  CONTINUE 

RETURN 

END 
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c**«.  **.•••*•.**•*•*.  *.***•»••»»• 

c 

C  SURF 
C 

C  THIS  ROUTINE  GENERATES  X,Y  LOCATIONS  OF  LEFT,  CENTER  AND 
C  RIGHT  ENDS  OF  GEOMETRY  PANELS  FOR  CYLINDER  WITH  OUTER  AND 
C  INNER  RADIUS 
C 

C  INPUT: 

C  NELEMS  -  NUMBER  OF  OUTSIDE  PANELS 

C  THE  SAME  NUMBER  OF  PANELS  IS  GENERATED  ON  INNER 

C  SURFACE 

C  RADOUT  -  OUTER  RADIUS 
C  RADIN  -  INNER  RADIUS 
C 

C  OUTPUT: 

C  ELEM (1,1), ELEM (1,2)  -  X,Y  OF  LEFT  END  OF  ELEMENT  J 
C  ELEM(I,3),ELEM(I,4)  -  X, Y  OF  RIGHT  END  OF  ELEMENT  J 
C 

SUBROUTINE  SURF 
implicit  real#8  (a-h,o-z) 

COMMON/GEOMET/NNODES , NELEMS , ORDER , EPS , NORM I , NORM J , JAC , NODEC , 
•NODPIN 1 , NODPIN  2 , NODE { 5 ) , ELEM( 360,4) 
integer  order 

real *8  NODE, NORM l , NORM J . JAC 


PI-3. 14159265 
READ ( 5, * )  RADOUT, RADIN 
1  FORMAT ( 2  F 1 0 . 4  ) 

GENERATE  FLAT  PANEL  GEOMETRY 

ANGLE-2. *PI /NELEMS 
PLENGO-2.*RADOUT*dTAN( ANGLE/2. ) 
PLENGI-2 . 'RADIN  *dTAN ( ANGLE/ 2 . ) 

CALCULATE  X,Y  LOCATIONS  FOR  OUTER  SURFACE 

DO  10  1-1, NELEMS 
ANG-U-1  )*ANGLE'PI 

RSIDE*SfiRT< (PLENGO/2. } *»2+RADOUT*,2 ) 
ELEM( 1,1 )-RSIDE*dCOS(ANG-( ANGLE/2. ) ) 
ELEM  (  1,2 )-RSIDE  *dSIN ( ANG- ( ANGLE/2 . ) ) 
ELEM( I , 3 ) -RSI DE*dC0S ( ANG+ ( ANGLE/2 . ) ) 
ELEM( 1,4 )«RSIDE  *dSlN ( ANG> ( ANGLE/2 .  ) ) 
10  CONTINUE 

CALCULATE  X,Y  LOCATIONS  FOR  INNER  SURFACE 
C 

DO  20  1-1, NELEMS 
ANG- ( 1-1 ) *ANGLE+PI 
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RSIDE«SQRT< (PLENGI/2. ) **2+RADIN**2  ) 
ELEM(NKLEMS+I, 1 )«RSIDE*dCOS ( ANG- ( ANGLE/2 . ) ) 
ELEM( NELEMS+1 , 2 )*RSIDE *dSIN { ANG- ( ANGLE/2 . ) > 
ELEM ( NE  LEMS+I , 3 ) “RSI DE*dCOS ( ANG+ ( ANGLE/2 . ) ) 
ELEMINELEMS+I , 4 )-RSIDE *dSIN ( ANG+ ( ANGLE/2. ) ) 
20  CONTINUE 

RETURN 
END 
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c 

C  VKLOC 
C 

C  THIS  SUBROUTINE  CALCULATES  THE  U  AND  V  COMPONENTS  OP  THE 
C  PERTURBATION  VELOCITY  (UVEL,WEL) 

C 

SUBROUTINE  VELOC 
implicit  real*8  (a-h,o-z) 

COMMON /FLAGS/IFLAG, PFLAG,  ITER 

COMMON /GEOMET/NNODES , NELEMS , ORDER , EPS , NORM I , NORM J , JAC , NODEC , 
•NODPIN 1,NODPIN2, NODE (5 ) ,ELEM( 360,4) 

COMMON/POSITION/XI ( 1080 ) ,YI (1080 ) 

COMMON/FLOW/H <240, 240 ) ,G( 240, 240) ,Q( 240), GHAT (240 ) 

COMMON /POTEN/VINF, ALPHA, VI ,VJ , PINF , RHOINF, GAMMA , M1NF2 , NODEP 
COMMON /VELOC  IT/UVEH 240 ) , WEL ( 240 ) , XFORC ( 480 ) , YFORC ( 480 ) , 

*XOLD( 480 ) , YOLD( 480 ) 

REAL *8  SHAP ( 5 ) , ETA( 12),WT(12),SHP(5) 

REAL*8  ETAS ( 8 ) , WTS ( 8 ) , SHP 1 ( S ) , SHP2 ( 5 ) , SHP3 ( 5  ) 

REAL *8  NODE, JAC, JAC 1, JAC 2, NORMI,NORMJ, GAMMA, MINF2 
INTEGER  ORDER, PFLAG 

REAL*8  TERMQX( 240 , 240 ) ,TERMQY (240, 240 ) , TERMPX (240,240), 

*  TERMPY ( 240, 240 ) 

C 

ETA  AND  WT  ARE  12TH  ORDER  GAUSSIAN  QUADRATRUE  VALUES  AND  WEIGHTS 

DATA  ETA  (  1  )  ,  ETA  (  2  )  ,  ETA  (  3  )  ,  ETA  (  4  )  ,  ETA  (  5  ) ,  ETA  ( 6  ) ,  ETA  (  7  ) , 

•ETA ( 8 ) , ETA ( 9 ) , ETA ( 10 ) , ETA ( 1 1 ) , ETA( 1 2 ) /-. 981 56063d0 , 904 1 1 726d0 , 
•- . 76990 2 6 7d0 ,-. 58 731795 dO,-. 36 7831 49d0 , 1252 334 IdO, . 1252 334 IdO, 
*. 36783 149d0, .58731 795d0,.76990267d0, .904 11 726d0, .981 56063d0/ 
DATA  WT ( 1 ) , WT ( 2 ) , WT ( 3 ) , WT ( 4 ) , WT ( 5 ) , WT ( 6  )  , WT ( 7 ) , WT ( 8 ) , WT ( 9 ) , 

*WT ( 10 ) ,WT (11) , WT (  12  )/. 04  7 1 753  dO , . 10693933d0, .  16007833d0, 
*.20316743d0, .233492 54d0, .24914705d0, .24914705d0, .23349254d0, 

• . 203 16743d0, . 1 600783 3d0, . 10693933d0, . 047 1 7534d0/ 

ETAS  AND  WTS  ARE  8TH  ORDER  GAUSSIAN  VALUES  AND  WEIGHTS  FOR 
1/R  SINGULAR  INTEGRALS 

DATA  ETAS ( 1 ) , ETAS ( 2 ) , ETAS ( 3 ) , ETAS ( 4 ) , ETAS ( 5 ) , ETAS ( 6 ) , ETAS ( 7 ) , 
•ETAS (8)/-.00324250dO,.05349077dO, . 1 7782733d0 , . 3 50 7 1 788d0 , 

• . 545819S2d0, .73342708d0, .88498305d0, .97745438d0/ 

DATA  WTS ( 1 ) , WTS ( 2 ) , WTS ( 3 ) ,  WTS ( 4 ) , WTS ( 5 ) , WTS ( 6 ) , WTS ( 7 ) , 

•WTS (8)/- 3. 93063681 dO, 1 . 7374063 IdO , .85763454d0, .53836077d0, 

*. 35975 759d0, .23726 769d0, . 14146228 dO, .05874 76 IdO/ 

DATA  PI/3. 14159265d0/ 

C 

C  INITIALIZE  MATRICES 
C 

DO  10  I-1.NNODES 
UVELII ) — 0 . 0 
WEL(I  )— 0 . 0 
DO  10  J“  1 , NNODES 
TERMPXf I , J )— 0 . 0 
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VELOC  2/6 


TERMPY ( I , J ) *0 . 0 
TERMQX( I , J  )«0 . 0 
TERMQY (I,J)“0.0 
10  CONTINUE 
C 

C  BEGIN  LOOPS  ON  NODE  I  AND  ELEMENT  J 

C  FOR  EACH  I  ROW,  PERFORM  INTEGRATIONS  OVER  EACH  ELEMENT  J 
C  AROUND  THE  STRUCTURE.  FILL  IN  ONE  ROW  AT  A  TIME 
C 

DO  1000  1*1 , NNODES 
DO  1000  J-1.NELEMS 

DETERMINE  IF  THE  INTEGRATION  SHOULD  BE  SINGULAR 
IE  IS  THE  NODE  I  ON  THE  ELEMENT  J? 

ISING  «  0  MEANS  NODE  IS  NOT  ON  ELEMENT 

ISING  -  1  MEANS  NODE  IS  FIRST  NODE  ON  ELEMENT  J 

ISING  «  2  MEANS  NODE  IS  SECOND  NODE  ON  ELEMENT  J 

ISING  -  3  ETC 

ISING-0 

DO  20  L»1, ORDER 

NODES-J ‘ORDER- ( ORDER-L ) 

IF( I . EQ. NODES )  ISING-L 
20  CONTINUE 

IF(ISING.EQ.O)  THEN 


THE  INTEGRATION  IS  NOT  SINGULAR 
PERFORM  1 2TH  ORDER  GAUSSIAN  QUADRATURE 


DO  50  K-1, 12 

CALCULATE  THE  X  AND  Y  COORDINATES  OF  THE  QUADRATURE  POINT 
AND  CALCULATE  DISTANCE  FROM  NODE  I  TO  QUAD  POINT 

ETAA—1 . 

ETAC- 1 . 

ETAQ-ETA  ( K  ) 

SHAP( 1  )" ( ETAQ-ETAC ) /{ ETAA-ETAC ) 

SHAP ( 2 )“ ( ETAQ-ETAA ) /( ETAC-ETAA ) 

X-SHAP(1 ) *ELEM ( J , 1)+SHAP(2) *ELEM( J , 3 ) 

Y-SHAP< 1 )*ELEM(J,2)+SHAP(2)*ELEM(J,4 ) 

RI-XI(I )-X 
RJ-YKI  )-Y 
R“SQRT( RI *RI+RJ  *RJ ) 

CALCULATE  DQ*/DX , DQ*/DY , DPH I */DX , DPHI */DY 
C 

CALL  JACOB I AN (J, ETA (K ) ) 

D“RI * ( NORMI ) +RJ  * ( NORMJ ) 

DPSTDX— RI/(  R*R) 

DPSTDY— RJ/(R*R) 

DQSTDX«( ( R*R* ( NORMI )-D*2 . *RI )/(R*R*R*R) ) 
DQSTDY-( ( R*R*( NORMJ )-D*2.*RJ)/(R*R*R*R) ) 

C 
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C  CALCULATE  DH/D  AND  DG/D  MATRIX  TERMS 
C  ROW  CORRESPONDS  TO  NODE  I 

C  COLUMN  CORRESPONDS  TO  GLOBAL  NUMBER  OF  ELEMENT  NODES 
C  I.E  NODE  1.2,3  ON  ELEMENT  1  ARE  1,2,3 

C  NODE  1,2,3  ON  ELEMENT  2  ARE  4,5,6 

C  NODE  1,2,3  ON  ELEMENT  3  ARE  7,8,9 

C 

CALL  SHAPE (ORDER, ETA ( K ) , SHAP ) 

DO  35  L-1, ORDER 

J J -J  *0  RDE  R  -  ( 0  RDE  R -L  ) 

TE RMQX ( I , J J )  -TE RMQX ( I , J J )  + 

*  SHAP ( L ) *DQSTDX*JAC*WT ( K ) /( 2 . *PI  ) 

TERMPX ( I , J J )  —TERMPX ( I , JJ  )♦ 

*  SHAP(L)*DPSTDX*JAC»WT(K)/{2.*P1 ) 

TERMQY ( I , JJ ) -TERMQY ( I , JJ )+ 

*  SHAP(L)*DQSTDY*JAC*WT(K)/(2.»PI  ) 

TERMPY ( I , J J ) -TERMPY ( I , J J ) + 

*  SHAP(L)*DPSTDY*JAC*WT(K)/(2.*PI ) 

35  CONTINUE 

50  CONTINUE 

ELSE 
C 

C*  •*••*****•******•*•••*•**••**•*• 

C  THE  INTEGRAL  IS  SINGULAR  ON  ELEMENT  J 
C 

C  THE  DPH I */D  TERMS  ARE  CALCULATED  SEPARATELY  FROM  THE  DQ*/  D  TERMS 
C  BECAUSE  DPHI*/D  IS  1/R  SINGULAR  AND  DQV D  IS  1/R»R  SINGULAR 
C 

C  CALCULATE  DQ»/D  TERMS  FOR  SINGULAR  ELEMENT 

C 

C 

C  BREAK  UP  INTERVAL  INTO  SINGULAR  AND  NON-SINGULAR  PART 
C  SINGULAR  PART  CENTERED  AROUND  SINGULARITY 
C  START-START  OF  SINGULAR  PORTION 
C  END-END  OF  SINGULAR  PORTION 
C 

SING-NODE (ISING) 

DELTA- ( 1. -ABS ( NODE ( ISING ) ) ) 

START-NODE (ISING) -DE  LTA 
END-NO DE( ISING  )  DELTA 
SUM-0 . 0 
ETAA— 1. 

ETA0 1 . 

C 

C  IF  START  OF  SINGULAR  INTEGRAL  IS  -1  AND  END  IS  1  THEN 
C  THE  ENTIRE  ELEMENT  IS  IN  THE  SINGULAR  INTERVAL 
C  SKIP  THE  NON-SINGULAR  CALCULATIONS 
C 

IF( START. EQ. -1 .dO .AND. END. Eg.  l.dO  )  GO  TO  300 
DO  250  K-1,12 
C 

C  START  OF  SINGULAR  INTERVAL  IS  -1,  END  IS  END 
C  THERFORE  NON-SINGULAR  PART  IS  FROM  END  TO  1 
C 
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IF< START. Eg. -1.  dO )  THEN 

ET-(ETA(K)+1. )*( 1.-ENDJ/2.+END 
CALL  SHAPE (ORDER, ET, SHP) 

ETAQ-ET 

SHAP( 1 )- ( ETAg-ETAC ) /( ETAA-ETAC ) 

SHAP ( 2 )- ( ETAQ-ETAA ) / ( ETAC-ETAA ) 

X«SHAP< 1 )*ELEM(J, 1 )+SHAP( 2 ) *ELEM( J , 3 ) 

Y-SHAP( 1 ) *ELEM( J, 2 )+SHAP(2 ) *ELEM( J,4 ) 

Rl-XI (I  )-X 
RJ-YI (I )-Y 
R-SgRT(RI*RI+RJ*RJ  ) 

CALL  JACOBIAN ( J , ET  ) 

DgSTDX-(NORMI )/( R*R ) 

DgSTDY- ( NO RMJ ) / ( R * R  ) 

DO  100  L-1, ORDER 

JJJ -J ‘ORDER- ( ORDER-L ) 

TERNgX ( I , J J J )=TERMgX ( 1 , J JJ ) + 

*  SHP ( L) ‘DQSTDX*JAC* . 5  * ( 1 . -END ) *WT ( K ) / ( 2 . *PI ) 
TERMgY ( I , JJJ ) -TERMgY ( I , JJJ )+ 

*  SHP ( L ) ‘DgSTDY  * JAC* . 5  * (  1 . -END ) *WT ( K )/ ( 2 . *PI ) 

100  CONTINUE 

ELSE 

C 

C  END  OF  SINGULAR  INTERVAL  IS  1 

C  THEREFORE  NON-SINGULAR  PART  IS  FROM  -1  TO  START 
C 

ET-(ETA(K)+1 .  )  *  (  START-*- 1  .  )/2.-1. 

CALL  SHAPE ( ORDER, ET, SHP) 

ETAg-ET 

SHAP ( 1 )  -  ( ETAQ-ETAC ) / ( ETAA-ETAC ) 

SHAP(2 )- (ETAQ-ETAA)/ (ETAC-ETAA) 

X«SHAP( 1 ) *ELEM ( J , 1 )+SHAP( 2 ) *ELEM( J , 3 ) 

Y-S HAP ( 1 ) ‘ELEM ( J , 2 ) +SHAP ( 2 ) ‘ELEM ( J , 4 ) 

RI-XKI  )-X 
RJ-YI  (I  ) -Y 
R-SgRT ( R1 *RI+RJ*RJ } 

CALL  JACOBIAN (J,ET) 

DgSTDX-(NORMI )/(R‘R) 

DgSTDY- (NO RMJ )/(R*R) 

DO  200  L-1, ORDER 

JJJ-J ‘ORDER- ( ORDER-L ) 

TERMgX ( I , J J J ) -TERMgX ( I , J J J ) + 

*  SHP(L)*DQSTDX*JAC‘.5*(START+1. )*WT(K)/(2.‘PI  ) 
TERMgY ( I , JJJ ) -TERMgY ( I , JJJ  )  + 

*  SHP (L ) ‘DgSTDY ‘JAC* . 5* ( START +  1  . ) *WT (K ) / ( 2 . *PI  ) 

200  CONTINUE 

END  IF 

250  CONTINUE 

C 

C  NOW  CALCULATE  SINGULAR  INTEGRAL 
C 

300  CONTINUE 

DO  400  K- 1,12 

X2X 1-ELEM( J , 3 )-ELEM< J , 1 ) 
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Y2Y 1 -ELEM( J, 4 ) -ELEM( J ,  2  ) 

PLENG-SQRT(X2X 1*X2X 1+Y2Y 1 *Y2Y 1 ) 

ET-. 5* ( END-START )*(ETA(K)+1. ) +START 
CALL  JACOB IAN  ( J ,  ET ) 

DQSTDX-NORMI *8. /( PLENG ‘PLENG* ( END-START)  ) 
DQSTDY -NORMJ*8 . /( PLENG ‘PLENG* ( END -START  )  ) 
CALL  SHAPE (ORDER, ET, SHP ) 

DO  400  L-1, ORDER 

JJJ-J*ORDER- (ORDER-L ) 

TEI&tQX  ( I ,  JJJ  )  -TERMQX  ( I ,  JJJ  )  + 

*  SHP ( L ) *DQSTDX*JAC*WT (K ) / ! ETA (K)“2*(2.*PI  )  ) 
TERMQY ( I , JJJ ) “TERMQY ( I , JJJ )  + 

*  SHP ( L ) ‘DQSTDY  *JAC*WT (K)/(ETA(K)**2*(2.*PI ) ) 

400  CONTINUE 

C 

C  CALCULATE  SUMMATION  OF  WT/ETA2 
DO  500  K-1,12 

SUM»SUM+WT ( K ) /ETA ( K ) *  *2 
500  CONTINUE 

CALCULATE  F ( 0 )  AND  ADD  -(SUM+2)F(0)  TO  TERMS 
C 

ET-. 5* ( END-START ) +STA  RT 
CALL  JACOB  IAN (J,ET) 

DQSTDX-NORM I *8 . / ( PLENG*PLENG * ( END-START ) ) 

DQSTDY -NORMJ  *8 . / ( PLENG  *PL£NG  * ( END-START  )  ) 

CALL  SHAPE (ORDER, ET, SHP ) 

DO  600  L-1, ORDER 

JJJ-J ‘ORDER- (ORDER-L ) 

TERMQX ( I, JJJ )— TERMQX ( I , JJJ ) - 

*  ( SUM+2 . ) *SHP(L )*DQSTDX*JAC/( 2 . *PI ) 

TERMQY ( I, JJJ ) -TERMQY (I, JJJ >- 

*  (SUM+2.  ) *SHP(L) ‘DQSTDY *JAC/( 2. *PI  ) 

600  CONTINUE 

C‘‘‘‘‘‘‘‘‘‘‘‘***‘*‘‘‘‘‘‘‘‘‘‘‘‘ 

J  CALCULATE  DPHI‘/D  TERMS  FOR  THE  SINGULAR  ELEMENT 
C 

X2X 1 -ELEM ( J , 3 ) -ELEM ( J, 1 ) 

Y  2Y 1 -ELEM ( J , 4 ) -ELEM ( J , 2 ) 

PLENG»SQRT(X2X 1 *X2X 1 +Y  2Y 1 *Y  2Y 1 ) 

PX-X2X 1  •  (  2 . /PLENG “2  )  * (  1./(2.*PI  ) ) 

PY -Y2Y  1  * (  2  .  /PLENG “2  )•(  1 ./(  2.  *PI  )  ) 

C 

C  PERFORM  8TH  ORDER  SUMMATION 
C 

DO  740  K-1,8 

ET-(-1 .-SING) ‘ETAS (KJ+SING 
CALL  JACOBIAN (J,ET) 

JAC 1 -JAC 

CALL  SHAPE (ORDER, ET , SHP  1  ) 

ET- ( 1. -SING) ‘ETAS (K)+SING 
CALL  JACOB IAN ( J , ET ) 

JAC2-JAC 

CALL  SHAPE (ORDER, ET,SHP2 ) 
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DO  730  L-1, ORDER 

JJJ -J ‘ORDER- (ORDER-L ) 

TERMPX ( I , JJ J )-TERMPX(I, JJJ)- 

*  SHP 1  (L  *PX*JAC  1  *WTS  (K  )  +SHP2 ( L )  *PX*JAC2*WTS  ( K  ) 
TERMPY ( I ,  JJJ ) -TERMPY ( I , JJJ )- 

*  SHP 1 ( L ) *PY • J AC  1 • WTS ( K.  1  +  SHP2  ( L  ) *PY  *JAC2 *WTS ( K ) 

730  CONTINUE 

740  CONTINUE 

CALCULATE  F ( S INGULARITY )  LOG  TERMS  AND  ADD  TO  THE  SUMMATION 

CALL  SHAPE (ORDER, SING, SHP3) 

CALL  JACOB  IAN ( J, SING ) 

DO  750  L-1, ORDER 

J JJ “J ‘ORDER- ( ORDER'L ) 

TERMPX ( I, JJJ ) -TERMPX ( I ,  JJJ  )- 

*  SHP3(L) *PX*JAC*dlog{ ABS( -1.-SING ) )+ 

*  SHP3(L)*PX*J  AC*  d log  (  ABS  (  1  .  -SING  )  ) 

TERMPY ( I, JJJ ) -TERMPY ( I, JJJ )- 

*  SHP3(L)*PY*JAC*dlog( ABS ( - 1 . -SING ) )+ 

*  SHP3 (L) *PY  *JAC*dlog( ABS (  1.-SING) ) 

750  CONTINUE 

END  IF 

1000  CONTINUE 

IF  PRINT  FLAG  -  1  PRINT  OLTT  TERM  MATRICES 

IF(PFLAG.E(2.  I.OR.PFLAG.EQ.  3}  CALL  VELOUT  ( TERMQX , 

*  TEWQY, TERMPX, TEfWPY) 

:  CALCULATE  U  AND  V  PERTURBATION  VELOCITIES 

DO  1200  I  —  1 , NNODES 

DO  1200  J-1, NNODES 

UVELQ  )— UVEL ( I  )  ♦  (  -TERMQX  ( I ,  J  >  *GHAT  ( J  ) 

*  +TERMPX(I,J)*Q(J) ) 

VVELU  )-WEL<I  )+(-TERMQY(I,J)*GHAT(J  ) 

*  +TERMPY ( I,J) *Q( J ) ) 

1200  CONTINUE 

DO  1250  1-1, NNODES 
UVEL(I)-UVEL(I)/.5 
WEL(  I  l-WELI  I  )/  .  5 
1250  CONTINUE 

do  1300  i-1,nnodes 
uvel( i  )«vi-uvel( i ) 
vvel ( i ) -vj-vvei ( i ) 

1300  continue 
RETURN 
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c 

C  VELOUT 
C 

C  THIS  SUBROUTINE  PRINTS  THE  INTERMEDIATE  PARTS  OF  THE 
C  VELOCITY  MATRICES  FOR  CHECKOUT 
C 

C** **************************************************************  * 

SUBROUTINE  VELOUT  ( TERMQX,  TEWQY ,  TERMPX,  TE5*iPY  ) 
implicit  real*8  (a-h,o-z) 

COMMON /GEOMET/NNODES , NELEMS , ORDER , EPS , NORM  I , NORMJ, JAC , NODEC , 
*NODPIN 1 , NODPIN2, NODE ( 5  ) , ELEM( 360,4) 
real*8  TERMQXf 240, 240 ) .TERMQY ( 240, 240 ) ,TERMPX( 240, 240  ) , 
•TERMPY (240,240) 
real *8  NORM  I , NORMJ , JAC , NODE 
INTEGER  ORDER 
C 
C 

WRITE (6, 1 ) 

I  FORMAT (//, 5X , 1 TERMQX  MATRIX*,/) 

DO  10  1-1, NNODES 

WRITE (6,5)  <TERMQX(I,J) ,J»1, NNODES) 

5  FORMAT ( /13(8(F10.4,5X) ,  / )  > 

10  CONTINUE 

WRITE( 6,11) 

II  FORMAT (//, SX, 'TERMQY  MATRIX*,/) 

DO  20  1-1, NNODES 

WRITE  ( 6-,  15)  (TERMQY  ( I ,  J ) , J*  1  ,  NNODES  ) 

15  FORMAT ( /13(8(F10.4,5X) ,  / ) ) 

20  CONTINUE 
WRITE( 6,21 ) 

21  FORMAT( //, 5X, ’TERMPX  MATRIX',/) 

DO  30  1-1, NNODES 

WRITE (6, 25)  (TERMPX< I, J ) , J-1, NNODES) 

25  FORMAT (/13(8(F10.4,5X) ,/) ) 

30  CONTINUE 
WRITE (6, 31 ) 

31  FORMAT (//, 5X , 'TERMPY  MATRIX’,/) 

DO  40  1-1, NNODES 

WRITE (6, 35)  (TERMPY ( I , J ) ,J*1, NNODES) 

35  FORMAT ( /13(8(F10.4,5X) ,/) ) 

40  CONTINUE 
RETURN 
END 


Appendix  C:  Background  Concepts 


Weighted  Residuals 

The  method  of  weighted  residuals  may  be  defined  in  the  following  manner.  Begin  with  the  problem 
definition; 

L(uq)  =  b  in  domain  Q  (Cl) 

with  boundary  conditions 

S(uq)  =  s  on  boundary  T  j  (C2) 

C(uq)  =  g  on  boundary  I~2 

where  u0  is  the  exact  solution  to  the  problem.  The  approximation  may  be  made  so  that 

u0  =  u  =  £  akNk  +  a0  (C3) 

k  >1 

where  arc  undetermined  coefficients  and  are  linearly  independent  functions.  Substituting  the 
approximation  into  Equations  Cl  and  C2: 

L(u)  -  b  *  0  =  R 

S(u)-s*0  =  R,  (C4) 

G(u)  -g  *  0  =  R2 

where  R,  Rj  and  R2  arc  the  residuals,  or  the  errors.  The  purpose  of  the  various  numerical  methods 

techniques  is  to  make  all  the  errors  as  small  as  possible  over  the  domain  and  boundary.  One  way  of  doing 
this  is  with  the  weighted  residual  method. 

Take  a  function  w,  such  that 

w  =  £  bkyk  (C5) 

k  » ! 
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The  error  can  be  distributed  throughout  domain  Q  by  multiplying  the  error  by  the  function  w  which  is 


called  a  weighting  function. 


R  wdfl  =  0 


(C  6) 


Because  the  b  terms  are  arbitrary,  the  key  is  to  choose  the  function  so  that 

Sa R  Vkd£l  =  o 


(C7) 


Shape  (Interpolation)  Functions 


The  value  of  a  function  u  may  be  defined  in  terms  of  its  value  at  known  points  (u^)  and  shape,  or 
interpolation,  functions  (1%)  such  that; 


uOl)  =  N i  (t|)  u ,  +  N2(T|)  u2  +  .  .  .  +  Nk(ri)uk 


(C8) 


k  is  the  order  of  the  function  such  as  1st  (constant),  2nd  (linear),  3rd  (parbolic),  4th,  etc.  The  points  at 
which  the  value  of  u  is  defined  are  called  the  nodal  points,  or  nodes. 

Lagrangian  polynomial  shape  functions  were  used  in  the  development  of  this  thesis  arc  may  be 
calculated  by: 


<M*1)  = 


(fi-fii).  ■  •  (fi-TUi)  ■  ■  (fi-fik) 

(fii-fii)  •  ■  •  (fi.-fi-i)  Oh-TUi)  •  •  • 


(C9) 


The  variable  i)  is  an  elemental  coordinate,  t)  varies  from  -1  at  the  left  end  to  +1  at  the  right  end  of  the 
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element.  The  nodes  are  therefore  distributed  throughout  the  element  with  values  from  -1  to  1.  For 
example,  a  third  order  function  may  have  nodes  defined  at  -1, 0,  +1. 


Jacohians 

A  Jacobian  allows  the  transformation  from  the  integral  about  the  boundary  (dr)  to  the  integral  along 
the  element  surface  (dri).  For  a  two-dimensional  problem  the  Jacobian  is  defined  as; 

111  -  V  (If + (If = I  (c,o> 

Therefore, 


dT  =  |J|  dn 


(Cll) 


In  order  to  calculate  the  Jacobian,  x  and  y  must  be  expressed  in  terms  of  the  variable  n>  For  a  flat  element 
where  (xj,  y})  are  the  coordinated  for  the  left-end  of  the  element  and  (x2,y2)  are  the  coordinates  for  the 
right-end,  x  and  y  may  be  expressed  by  shape  functions; 

x  =  Nj  Xj  +  N2  x2  y  =  Nj  yj  +  N2  y2  (C12) 


and 


ax 

an  ap 


+ 


ad>2 

dn 


X2 


dy  dOj 
2  ~ —  yi 

dn  dn 


(Cl  3) 
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